%% cia_3e.m
%
%  Basic cash-in-advance model from Chapter 3, 
%  Monetary theory and Policy, Third edition,
%  Carl E. Walsh
%  MIT Press, copyright
%
%  See cia_approx.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 = 7; % z, u, y, n, k, m, i
n2 = 2; % pi lambda

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;      % 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

%% 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;    % 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)^(-p.alpha/(1-p.alpha)); % = 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.alpha p.delta p.beta p.PHI p.eta p.a p.b  ]
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), n(t), k(t), m(t), i(t)
%  E_t(pi(t+1)), E_t(lambda(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,4) = p.alpha - 1;
A2(3,1) = 1;
A2(3,5) = p.alpha;

% labor supply
A1(4,3) = -1;
A1(4,4) = (1+p.kappa);
A2(4,9) = 1;

% capital accumulation/resource constraint)
A1(5,3) = -p.YK_bar;
A1(5,5) = 1;
A1(5,7) = p.CK_bar;
A2(5,5) = 1 - p.delta;

% Marginal utility of consumption and income(lambda)
A1(6,6) = -(p.i_bar/p.I_bar);
A1(6,7) = -p.PHI;
A2(6,9) = 1;

% money stock
A1(7,7) = 1;
A2(7,2) = 1;
A2(7,7) = 1;
A2(7,8) = -1;

% Euler condition
A1(8,6) = 1;
A1(8,8) = -1;
A1(8,9) = 1;
A2(8,9) = 1;

% real return
A1(9,5) = -p.alpha*p.YK_bar*(1-p.alpha)*p.kappa/(p.alpha + p.kappa);
A1(9,6) = -1; 
A1(9,8) = 1;
A1(9,9) = p.alpha*p.YK_bar*(1-p.alpha)*p.kappa/(p.alpha + p.kappa);
A2(9,1) = -((1+p.kappa)/(p.alpha + p.kappa))*p.rho_z;

A = 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] = olrp_tr(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), n(t), k(t), i(t) m(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)]'

%disp('Standard deviation of output')
%sdy = sqrt(Sx1(3,3))
%disp('Standard deviations relative to output')
%sdv = sqrt(diag(VarCov));
%sdv./sdy
%% 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('shock, output, employment, capital, nominal rate, and inflation')
        disp('Hit any key to continue')
        pause
        resp_z = [x(st:Tbig-1,1) x(st+1:Tbig,3:6) x(st+1:Tbig,8)  ]; 
        resp_z
        figure (1)
        title('Response to a productivity shock')
        subplot(2,1,1)
        plot(tid,resp_z(:,2),tid,resp_z(:,4),tid,resp_z(:,3),tid,tzero,'k')
        legend('output','capital','employment')
        subplot(2,1,2)
        plot(tid,resp_z(:,5),tid,resp_z(:,6),tid,tzero,'k')
        legend('nominal interest rate','inflation')
    else
        disp('Responses to a money growth shock')
        disp('shock, output, employment, capital, nominal rate, and inflation')
        disp('Hit any key to continue')
        pause
        resp_v = [x(st:Tbig-1,1) x(st+1:Tbig,3:6) x(st+1:Tbig,8)  ]; 
        resp_v
        figure (2)
        title('Response to a money growth rate shock')
        subplot(2,1,1)
        plot(tid,resp_v(:,2),tid,resp_v(:,3),':',tid,tzero,'k')
        legend('output','employment')
        subplot(2,1,2)
        plot(tid,resp_v(:,5),tid,resp_v(:,6),':',tid,tzero,'k')
        legend('nominal interest rate','inflation')
    end
end; % end of loop over shocks (sk)
end; % end of loop over imp
resp3_v = resp_v;