% BASIC_RBC_ps2.M:
% Used for problem set 2, winter, 2005
%
% choice of utility of leisure function
% for (n^(1+eta))/(1+eta), utility = 1;
% for ln(1-n), utility = 2;
utility = 2;

% Setting parameters:

Z_bar     = 1;    % Normalization
alpha     = .66;  % Labor's share
delta     = .025; % Depreciation rate for capital
R_bar     = 1.01; % One percent real interest per quarter
N_bar     = 0.33; % fraction of time in market work
sig       = 1.0;  % constant of relative risk aversion = 1/(coeff. of intertemporal substitution)
eta       = 1.0; 
rho       = .2;  % autocorrelation of technology shock
sigma_eps = .006; % Standard deviation of technology shock.  Units: Percent.

% Calculating the steady state:

beta    = 1.0/R_bar;  % Discount factor beta
YK_bar  = (R_bar -1 + delta)/(1-alpha);
CK_bar  = YK_bar - delta;
CY_bar  = CK_bar/YK_bar;

if utility == 1;
    le = eta;
else
    le = N_bar/(1-N_bar);
end
% Declaring the matrices. 
    % The first thing to do is name the variables
    % Elements in VARNAMES must be of same size so if
    % some names are shorter than others, fill in with blanks. 

VARNAMES = ['capital    ',
            'output     ',
            'consumption',
            'investment ',
            'employment ',
            'return     ',
            'technology '];

% Translating into coefficient matrices.  
% The loglinearized equations are, conveniently ordered:
% 1) 0 = K(t+1) - (1 - delta) k(t) - delta I(t)
% 2) 0 = Y(t) - (1-a) K(t) - a N(t) - z(t)
% 3) 0 = Y(t) - sig C(t) - (1 + le) N(t)
% 4) 0 = Y(t) - (C/Y) C(t) - delta*(K/Y) I(t) 
% 5) 0 = - beta*(1-alpha)*Y/K [Y(t) - K(t)] + R(t) 
% 6) 0 = EC(t+1) - C(t) - (1/sig)ER(t+1)
% 7) z(t+1) = psi z(t) + epsilon(t+1)
% CHECK: 7 equations, 7 variables.
%
% Endogenous state variables "x(t)": k(t)
% Endogenous other variables "y(t)": y(t), c(t), i(t), n(t) r(t)
% Exogenous state variables  "z(t)": z(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,

% DETERMINISTIC EQUATIONS:

% for k(t+):
AA = [ 1
   0
   0
   0
   0   ];

% for k(t):
BB = [ -(1-delta)
    -(1 - alpha)
   0
   0
   beta*(1-alpha)*YK_bar  ];

% For [ y(t)  c(t), i(t)  n(t) r(t) ]
CC = [     0,   0,  -delta,  0,    0   % Equ. 1)
           1,   0,    0  -alpha,   0   % Equ. 2)
           1,  -sig,   0,  -(1+le), 0   % Equ. 3
           1   -CY_bar, -delta/YK_bar, 0  0   % Equ. 4
    -beta*(1-alpha)*YK_bar 0,  0,   0,  1]; % Equ. 5)

% For z(t):
DD = [ 0
   -1
   0
   0
   0 ];

% EXPECTATIONAL EQUATIONS:

% For k(t+1)
FF = [ 0 ];

% For k(t)
GG = [ 0 ];

% For k(t-1)
HH = [ 0 ];

% For [ y(t+1), c(t+1), i(t+1), n(t+1), r(t+1) ]
JJ = [ 0,  1,  0,  0,  -1/sig     ];

% For [ c(t), r(t), y(t) ]
KK = [ 0, -1,  0,   0,    0  ];

% For z(t+1)
LL = [ 0 ];

% For z(t)
MM = [ 0 ];

% AUTOREGRESSIVE MATRIX FOR z(t)

NN = [rho];

Sigma = [ sigma_eps^2  ];


% Don't chance the following--the is how the program knows how many of
% %each type of variable is in your model
[l_equ,m_states] = size(AA);
[l_equ,n_endog ] = size(CC);
[l_equ,k_exog  ] = size(DD);

% Setting the options:
HORIZON = 48; %how far out should the inpulse response be calculated
PERIOD     = 4;  % number of periods per year, i.e. 12 for monthly, 4 for quarterly
GNP_INDEX  = 2; % Index of output among the variables selected for HP filter
%  IMP_SELECT is a vector containing the indices of the variables to be plotted
IMP_SELECT = [2 3 5 7]; % plots vars. 2, 3, 5, 7 (y,c,n,e)
%IMP_SELECT = 1:(m_states+n_endog+k_exog); % plots all variables   
HP_SELECT  = 1:(m_states+n_endog+k_exog); % Selecting the variables for the HP Filter calcs.
DO_MOMENTS = 1;
DISPLAY_IMMEDIATELY = 1;
DISPLAY_AT_THE_END = 1;
DO_SIMUL = 0;
SIM_MODE = 1;          % Set to = 1, if you just want one simulated time series.
%                      % Set to = 2, if you want to do simulation-based calculation
%                      % of moments. 
SIM_LENGTH  = 100;     % Length of time series to be simulated. 
SIM_RANDOM_START = 1; 
%                      % = 1, if you want a random start, = 0, if you want predetermined start
SIM_GRAPH = 1;

% do_it calls the programs that actually solve the model
do_it;