clear all

% miu_3e.m: used to solve MIU model in Chapter 2 of C.E.Walsh, "Monetary 
% Theory and Policy," second edition."
% This model is solved using Harald Uhlig's Toolkit.
% This program is the same as MIU_2e.m from the second edition except for 
% some additonal notes.
disp(' Money in the utility function model,');
disp(' used in Chapter 2 of Monetary Theory and Policy 2nd ed.');
disp(' June 2001');


%% 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 =  2.56;      % interest elasticity
p.a = 0.95;       %
p.THETA = 1.0125; % rate of money growth (gross)

%% 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.0125; % rate of money growth (gross)
p.sigma_m = .95;% standard deviation of money growth rate shock. Units: Percent
p.rho_u = 0.67;    % autocorrelation of money shock
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 /(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]


p

disp('hit any key to continue')
pause

%Declaring the matries.

VARNAMES = ['capital        ',
            'money balances ',
            'output         ',
            'consumption    ',
            'employment     ',
            'lambda         ',
            'nominal rate   ',
            'real rate      ',
            'inflation      ',
            'investment     ',
            'technology     ',
            'money growth   '];
    
%Translating into coefficient matries.

%Check: 12 variables.
%Endogenous state variables "x(t)": k(t), m(t)
%Endonenous other variables "y(t)": y(t), c(t), n(t), lambda(t), i(t), r(t), inf(t), x(t).
%Exogenous state variables "z(t)": z(t), g(t).
%Switch to that notation. Find matrices for format
%0 = AA x(t) + BB x(t-1) + CC y(t) + DD z(t)
% 0 = E_t [FF x(t+1) + GG x(t) + HH x(t-1) + JJ y(t+1) ++ KK y(t) = LL z(t+1) + MM z(t)]
% z(t+1) = NN z(t) + epsilon (t+1) with E_t [ epsilon(t+1) ] = 0,
%
% Solution takes the form
% x(t) = PP*x(t-1) + SS*z(t)
% y(t) = RR*x(t-1) + QQ*z(t)

%equation order
%1). resource constraint
%2). production function
%3). money demand equation
%4). evolution of m
%5). labor market
%6). euler equation 1 
%7). euler equation 2
%8). investment evolvement equation
%8). fisher equation
%9). ex-post real interest

% for k(t), m(t):

AA = [1, 0
    0, 0
    0, 1
    0, 1
    0, 0
    0, -p.OMEGA2
    -p.alpha *(1-p.alpha)* p.YK_bar*(p.eta*p.NL_bar), 0
    1, 0];

%for k(t-1), m(t-1):
BB = [ (p.delta -1),  0
    p.alpha, 0
    0, 0
    0, -1
    0, 0
    0, 0
    0, 0
    (p.delta-1), 0];

%order: y(t), c(t), n(t), lambda(t), i(t), r(t), inf(t), x(t)
CC = zeros(8, 8);
CC(1, 1) = -p.YK_bar;
CC(1, 2) = p.CK_bar;
CC(2, 1) = -1;
CC(2, 3) = 1-p.alpha;
CC(3, 2) = -1;
CC(3, 5) = (1/p.b)*(1/p.i_bar);
CC(4, 7) = 1;
CC(5, 1) = 1;
CC(5, 3) = -(1 + p.eta*p.NL_bar);
CC(5, 4) = 1;
CC(6, 2) = -p.OMEGA1;
CC(6, 4) = 1;
CC(7, 4) = p.alpha*(1-p.alpha)*p.YK_bar;
CC(7, 6) = -p.theta;
CC(8, 8) = -p.delta;

DD = zeros(8, 2);
DD(2, 1) = 1;
DD(4, 2) = -1;
DD(7, 1) = p.alpha*p.YK_bar*(1+p.eta*p.NL_bar)*p.rho_z;

FF = zeros(2, 2);

GG = zeros(2, 2);

HH = zeros(2, 2);

JJ = zeros(2, 8);
JJ(1, 7) = 1;
JJ(2, 4) =1;

KK = [0, 0, 0, 0, -1, 1, 0, 0
    0, 0, 0, -1, 0, 1, 0, 0];

LL = [0, 0
    0, 0];

MM = [0, 0
    0, 0];

NN = [p.rho_z, 0
    p.phi, p.rho_u];

Sigma = [p.sigma_z^2 0
    0 p.sigma_v^2];

%setting the options:
[l_equ, m_states] = size(AA);
[l_equ, n_endog] = size(CC);
[l_equ, k_exog] = size(DD);

HORIZON = 32; %how far out should the inpulse response be calculated
PERIOD = 4; %number of periods per year, i.e. 12 for monthly, 4 for quarterly
DO_PLOTS = 0; %if impulse response plots should be made, = 0, if not
% IMP_SELECT is a vector containing the indices of the variables to be plotted
IMP_SELECT = [ 3, 5, 6, 7]; % for money growth shock
%HP_SELECT = 1: (m_states + n_endog + k_exog); %selecting the variables for the HP filter calcs.
%GNP_INDEX = 1; %Index of output among the variables selected for HP filter
DO_MOMENTS = 1;
DISPLAY_AT_THE_END = 0;
%Start the calculations;
do_it;
if DO_DISP3,
    disp('Standard deviations, HP-filtered series relative to output, cov with output:');
    VARNAMES
    disp([sqrt(varvec_fil)./sqrt(varvec_fil(3)) autcor_fil(1:12,6)]);
end;