%% miu_model_3e.m
%
%  Basic Money-in-the-utility function model from Chapter 2, 
%  Monetary theory and Policy, Third edition,
%  Carl E. Walsh
%  MIT Press, copyright
%
%  See miu_dynamics.tex for more details
% 
%  If sol == 'SHU', solution method basic on Toolkit of Paul Soderlind
%  If sol == 'LRP', solution method basic on linear regular approach
% 
%  This program sets up the basic model in the form
%  A1*E_t(Z_t+1) = A2 * Z_t + V_t+1
%
%     E_t(Z_t+1) = inv(A1)*A2 * Z_t + inv(A1)*V_t+1
%
%     E_t(Z_t+1) = A * Z_t + PSI_t+1
%
%  Solution takes the form
%
%     Z1(t+1) = M*Z1(t) + shocks
%
%     Z2(t) = C*Z(t)
%
%  where Z1 is the vector of predetermined variables and Z2 is the vector
%  of forward-looking variables.
%
%  Actual solution calculated by basicsol.m if sol == 'SHU'
%  [M,C] = basicsol(A,n1,n2,SigmaXX,x10,cutoff);
%  or if sol == 'LRP' by
%  [M, G] = lrp(A,p.beta,n1,n2)
%
%  basicsol.m (or lrp.m) must be on path
%
%  n1 = number or predetermined variables
%  n2 = number of nonpredetermined (forward looking) variables
%  SigmaXX = (n1+n2) x (n1+n2) covariance matrix Cov[PSI(t)] of shocks.
%  x10 is an n1x1 vector, initial values for predetermined variables
%  cutoff = scalar, max modulus of backward looking variables (1.0)
% 
%% Housekeeping
clear all
dispm = 0; % if = 1, matrices A1, A2, and A are displayed, if = 0, they are not displayed
sol = 'SHU'; % or 'SHU' or 'LRP'
imp = 1;  % set to 1 to calculate impulse resonses, 0 otherwise
%% Setting up model
n1 = 10; % z, u, y, c, k, x, n, r, m, i
n2 = 2; % lambda, pi

n = n1 + n2;
A1 = zeros(n,n);
A2 = zeros(n,n);

%% Parameters

% "Technology" parameters
p.N_bar = 1/3;    % steady state employment is a third of total time endowment
p.alpha = .36;    % Capital share
p.delta = .019;   % Depreciation rate for capital

%"Taste" parameters
p.beta = .989;    % discount parameters
p.PHI = 2;        % coefficient of relative risk aversion
p.eta = 1;        % leisure parameter

% Next two parameters are used for the MIU model but not for the CIA model
p.b = 3; % 5; %40; %2.56;      % interest elasticity
p.a = 0.99;       %

%% Exogenous stochastic processes

p.rho_z = .95;      % autocorrelation of techonology shock
p.sigma_z = .34; % standard deviation of techonology shock. Units: percent
    % 0.34 for PHI=2 and .48 for PHI = 8.

%% Exogenous policy specification
p.THETA = 1.01; % rate of money growth (gross)
p.rho_u = 0.5;    % autocorrelation of money shock
p.sigma_m = (.9/((1-(p.rho_u)^2)))^.5;% standard deviation of money growth rate shock. Units: Percent
p.phi = 0; %.15;    % correlation of techonology shock and money shock

% Innovation variance of money growth shock set so that s.d. of money
% growth matches p.sigma_m.
p.sigma_v = ((1-p.rho_u^2)*(p.sigma_m)^2 - ((p.phi^2)/(1-p.rho_z^2))*(p.sigma_z^2))^.5;

%% Calculating the steady state:
p.R_bar = 1/p.beta;                         % one percent real interest per quarter
p.YK_bar = (p.R_bar - 1 + p.delta)/p.alpha; % = Y_bar/K_bar
p.CK_bar = p.YK_bar - p.delta;              % = C_bar/K_bar
p.NK_bar = (p.YK_bar)^ (1/(1-p.alpha));     % = N_bar/K_bar
p.inf_bar = p.THETA - 1;                    % = inflation rate
p.NL_bar = p.N_bar/(1.0 - p.N_bar);         % = N_bar/L_bar
p.YN_bar = p.YK_bar /(1.0 - p.N_bar);       % = Y_bar/N_bar
p.I_bar = p.R_bar*p.THETA;                  % = (gross) nominal interest rate
p.i_bar = p.I_bar - 1;                      % = nominal interest rate

% Following used in MIU model but not in CIA model
p.mK_bar = (p.CK_bar) * (p.a*p.i_bar/((1-p.a) * (1+p.i_bar)))^ (-1/p.b); % = M_bar/K_bar
p.mC_bar = p.mK_bar /p.CK_bar;   % = m_bar/C_bar
p.H = (p.a*(p.CK_bar)^(1-p.b)) + (1-p.a)*(p.mK_bar^(1-p.b));
p.gamma = (p.a*(p.CK_bar)^(1-p.b))/p.H;

% Parameters defined in text
p.kappa = p.eta*p.N_bar/(1-p.N_bar); 
p.xi = p.alpha*(p.YK_bar)/(p.alpha + p.kappa + (1-p.alpha)*p.alpha*(p.YK_bar));
p.OMEGA1 = p.b*(p.gamma-1) - p.gamma*p.PHI;
p.OMEGA2 = (p.b - p.PHI)*(1-p.gamma);
p.theta = p.alpha + p.eta*p.NL_bar + (1-p.alpha)*p.alpha*p.YK_bar;

disp('Table 2.1: baseline parameter values')
disp('p.alpha p.delta p.beta p.PHI p.eta p.a p.b')
[p.a p.b p.alpha p.delta p.beta p.PHI p.eta ]
disp(' p.THETA p.rho p.sigma_z p.sigma_m')

[ p.THETA p.rho_z p.sigma_z p.sigma_m]
disp('Table 2.2: Steady-state values')
disp('[R_bar YK_bar CK_bar p.mK_bar NK_bar]')
[p.R_bar p.YK_bar p.CK_bar p.mK_bar p.NK_bar]
disp('hit any key to continue')
pause

%% Model matrices
% equations for z(t+1), u(t+1), y(t), c(t), k(t), x(t) n(t), r(t), m(t), i(t)
% E_t(lambda(t+1)), E_t(pi(t+1))

% exogenous disturbances
A1(1,1) = 1;
A2(1,1) = p.rho_z;
A1(2,2) = 1;
A2(2,1) = p.phi;
A2(2,2) = p.rho_u;

% production function
A1(3,3) = 1;
A1(3,7) = p.alpha - 1;
A2(3,1) = 1;
A2(3,5) = p.alpha;

% resource constraint (goods market clearing)
A1(4,3) = p.YK_bar;
A1(4,4) = -p.CK_bar;
A1(4,6) = -p.delta;

% capital accumulation
A1(5,5) = 1;
A1(5,6) = -p.delta;
A2(5,5) = 1 - p.delta;

% labor supply
A1(6,3) = -1;
A1(6,7) = (1+p.kappa);
A2(6,11) = 1;

% Marginal utility of consumption (lambda)
A1(7,4) = p.OMEGA1;
A1(7,9) = p.OMEGA2;
A2(7,11) = 1;

% real return
A1(8,5) = (1-p.alpha)*p.alpha*p.YK_bar*p.eta*p.NL_bar;
A1(8,8) = p.theta; 
A2(8,1) = p.alpha*p.YK_bar*(1 + p.eta*p.NL_bar)*p.rho_z;
A2(8,11) = (1-p.alpha)*p.alpha*p.YK_bar;

% money stock
A1(9,9) = 1;
A2(9,2) = 1;
A2(9,9) = 1;
A2(9,12) = -1;

% money demand
A1(10,4) = -1;
A1(10,9) = 1;
A1(10,10) = (1/p.b)*(1/p.i_bar);

% Euler condition
A1(11,8) = 1;
A1(11,11) = 1;
A2(11,11) = 1;

% Fisher equation
A1(12,8) = 1;
A1(12,10) = -1;
A1(12,12) = 1;

A = inv(A1)*A2

% Displaying the matrices that define the model
if dispm == 1;
disp('A1 matrix')
A1
disp('Hit any key to continue')
pause

disp('A2 matrix')
A2
disp('Hit any key to continue')
pause

disp('A, where [X(t+1) Ex(t+1)] = A*[X(t) x(t)] + [e(t+1) 0]')
A
end

% Variance-covariance matrix of shocks
SigmaXX = zeros(n1,n1);
SigmaXX(1,1) = p.sigma_z^2;
SigmaXX(2,2) = p.sigma_v^2;

Tbig = 24;  % length of impulse response
st=1;

%% Solving the model
x10 = zeros(n1,1);                  %%initial state vector*/

% check eigenvalues
rtck = 0;
[Xi_eigvec,Xi_eigval] = eig(A);
[Xi_sortabs,Xi_sortindex] = sort(abs(diag(Xi_eigval)));
if Xi_sortabs(n1+1) < 1
    rtck = 1;
    disp('Not enough unstable roots')
else
    if Xi_sortabs(n1) > 1
        rtck = 1;
        disp('Not enough stable roots')
    end
end

%% Solving the model
if rtck == 0
    if sol == 'SHU'; % solution method via Shur decomposition
        tic
        [M,C] = basicsol(A,n1,n2,SigmaXX,x10,1.0);
        toc
        disp('Solution obtaining using basicsol.m')
        disp('M, where X(t+1) = M*X(t) + e(t+1) for predetermined variables')
        M
        disp('C, where x(t) = C*X(t) for forward-looking variables')
        C
    else
        if sol == 'LRP' % solution method using linear regular approach
            tic
            [M, C] = lrp(A,p.beta,n1,n2)
            toc
            disp('Solution obtained using olrp_tr.m')
            disp('M, where X(t+1) = M*X(t) + e(t+1) for predetermined variables')
            M
            disp('C, where x(t) = C*X(t) for forward-looking variables')
            C
        end
    end
    disp('hit any key to continue')
    pause
end; % end of rtck loop

%% Calculating moments
% SigmaZ1 = M'*SigmaZ1*M + SigmaXX;
i = 1;
Sx10 = zeros(n1,n1);        %%iteraing on matrix*/
Sx1 =Sx10;
Vdiff = 100;
while Vdiff > 1E-6;
    Sx1 = M*Sx10*M' + SigmaXX;
    Vdiff = max(  max( abs(Sx10 - Sx1) )  );
    Sx10 = Sx1;
    i = i + 1;
end;

Sx2 = C*Sx1*C';  % Variance-cov. matrix of forward-looking variables
disp('variance-covariance matric')
VarCov = [Sx1 Sx1*C'; C*Sx1 C*Sx1*C'] ;

% Alternative methods for calculating variance-covariances of n1 variables
%SSV = vcv(SigmaXX,M,eye(n1)); 
%SV1=Sigkk(M,SigmaXX);

disp('Standard deviation of output')
sdy = sqrt(Sx1(3,3))
% equations for z(t+1), u(t+1), y(t), c(t), k(t), x(t) n(t), r(t), m(t), i(t)

disp('Standard deviations, s.d. relative to output, corr with output')
disp('Order of variables')
disp('y(t), n(t), c(t), k(t), x(t) r(t), m(t), i(t), lambda(t), pi(t)')
sdv = sqrt(diag(VarCov));
covarwithy = VarCov(3,:)';
corrwithy = covarwithy(:,1)./(sdv(:,1)*sdv(3,1));
[sdv(3,1) sdv(7,1) sdv(4,1) sdv(6,1) sdv(8:10,1)' sdv(12,1);
    sdv(3,1)/sdv(3,1) sdv(7,1)/sdv(3,1) sdv(4,1)/sdv(3,1) sdv(6,1)./sdv(3,1) sdv(8:10,1)'./sdv(3,1) sdv(12,1)/sdv(3,1)
    corrwithy(3,1) corrwithy(7,1) corrwithy(4,1) corrwithy(6,1) corrwithy(8:10,1)' corrwithy(12,1)]'

%% impulse responses
if imp == 1;
    disp('Impulse response functions')
else
        disp('Impulse reponses not chosen')
end
%% Constructing impulse responses
if imp == 1;
for sk = 1:2;  % loops over both shocks to construct impulse response functions
    if sk == 1;
        Shock0 = [1; zeros(n1-1,1)];   %%Impulse response wrt technology innovation*/
    else;
        Shock0 = [0; 1; zeros(n1-2,1)];   %%Impulse response wrt money innovation*/
    end
    
% Borrowed from Paul Soderlind's VarImpResp.m
x1 = zeros(Tbig,n1);   %%Simulate impulse response function*/
x1(:,:) = NaN;
x1M_t_1 = zeros(n1,1);           %%starting vector*/
t = 1;
  while t <= Tbig;
    if t == 1;
      Shock = Shock0;
    else;
      Shock = zeros(n1,1);
    end;
    x1M = M*x1M_t_1 + Shock;
    x1(t,:) = x1M';
    x1M_t_1 = x1M;
    t = t + 1;
  end;

%x1 = VarImpResp( M,Shock0,Tbig );  %VAR of x1(t)*/
    x2 = x1*C';                        %x2(t)*/
    x  = [x1 x2];                      %x(t) = [x1(t),x2(t)]*/
    
    tid = linspace(0,Tbig-1,Tbig-1)';
    tzero = tid - tid;
    % save results of impulse responses
    if sk == 1;
        disp('Responses to a productivity shock')
        disp('Hit any key to continue')
        pause
        resp_z = [x(st:Tbig-1,1) x(st+1:Tbig,3:8) x(st+1:Tbig,10) x(st:Tbig-1,12)  ]; 
        disp('z, output, cons., capital, investment, employment, real rate, nominal rate, and inflation')
        resp_z
        figure (1)
        subplot(2,1,1)
        plot(tid,resp_z(:,2),tid,resp_z(:,4),tid,resp_z(:,6),tid,tzero,'k')
        legend('output','capital','employment')        
        title('Responses to a productivity shock')
        subplot(2,1,2)
        plot(tid,resp_z(:,7),tid,resp_z(:,8),tid,resp_z(:,9),tid,tzero,'k')
        legend('real interest rate','nominal interest rate','inflation')
    else
        disp('Responses to a money growth shock')
        disp('Hit any key to continue')
        pause
        resp_v =  [x(st:Tbig-1,2) x(st+1:Tbig,3:8) x(st+1:Tbig,10) x(st:Tbig-1,12)  ]; 
        disp('u shk, output, cons., capital, investment, employment, real rate, nominal rate, and inflation')        
        resp_v
        figure (2)
        subplot(2,1,1)
        plot(tid,resp_v(:,2),tid,resp_v(:,6),':',tid,tzero,'k')
        legend('output','employment')        
        title('Responses to a money growth shock')
        subplot(2,1,2)
        plot(tid,resp_v(:,7),tid,resp_v(:,8),tid,resp_v(:,9),':',tid,tzero,'k')
        legend('real int. rate','nominal interest rate','inflation')
    end
end; % end of loop over shocks (sk)
end; % end of loop over imp
