%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% created 20 February 2017: 
% To solve the model with fiscal policy
% Wage and rental rate forecasts
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main file: calibrating shock standard deviation
% -  generate shock sequence, define all parameters, steady state 
%    values here and define them as global variables (unlike in EP)
% -  load the RE coefficients and define them as global (more efficient than EP)
% -  set up the minimization (call the optimization procedure: Sim's or fminunc Matlab)
%    *  input: initial values including ?_{A} (possibly Hessian matrix)
%    *  output: optimal ?_{A}
%
% Using Dynare solution: must ensure the parameters are the same
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
close all;
clc;

Font_size = 11;  
Line_width = 1.5; 
load coeff_k;
load coeff_u; 
display('start')

global delta alpha gamma eps_H tauk tauH tauC gy_ratio sigma beta_tilda R_bar beta ...
     theta ik_ratio ky_ratio iy_ratio by_ratio ay_ratio cy_ratio psi lam1 lam2 ...
     sigma_x   R_tilda epsc epsw kai gamma1 gamma2 ... 
     rk wage output rt hours cons invst caput cap gamma_x  gt at bt...
     epsZ  shockmatx   gain_control...
     coeff_mat_k_AL_ini coeff_mat_b_AL_ini   coeff_mat_R_AL_ini   coeff_mat_w_AL_ini...
      coeff_vec_k_RE coeff_vec_s_RE     ...
    phimatrix_RE omegamatrix_RE cvec_RE muvec_RE vecsigma_RE Mom_mat_ini ini_wR mom_mat_RE coeff_vec_g_RE

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

delta = 0.069/4;  %% depreciation rate

alpha = 0.41;   %0.34; %% capital share

gamma = exp(0.02/4);  %exp(0.004);%  % gamma_bar in the paper

eps_H = 1/10;%
%1/10000000000000;  % inverse of frisch elasticity

%tauk = 0.36;
tauk = 0.35;

%tauH = 0.4;
tauH = 0.46;

tauC = 0.18;  % 0.17
gy_ratio = 0.27; 
 
sigma = 2; %% parameter in the utility function

beta_tilda = 0.995; 

R_bar = gamma/beta_tilda; 

beta = gamma^sigma/R_bar;

theta = (R_bar-1)/(1-tauk)/delta+1;

ik_ratio = 1-(1-delta)/gamma; %% investment to capital ratio

ky_ratio = alpha*gamma/((R_bar-1)/(1-tauk)+delta);

iy_ratio = ky_ratio*ik_ratio;

%by_ratio = (-alpha*tauk-(1-alpha)*tauH+gy_ratio+delta/gamma*ky_ratio*tauk)/(1-R_bar/gamma);
by_ratio =  0.6;%0.75  0.728; 

ay_ratio = by_ratio+ky_ratio;

cy_ratio = 1-iy_ratio-gy_ratio;

psi = (1-alpha)/cy_ratio*(1-tauH)/(1+tauC);

lam1 =  0;

lam2 = -0.02;

sigma_x = 0.003;




R_tilda= beta*theta*delta*(1-tauk)/gamma^sigma/beta_tilda; 


%%% newly added parameters for models with wage forecasts


epsc = (1+tauC)*cy_ratio/ay_ratio + (1-alpha)/ay_ratio*(1-tauH)/(eps_H-(sigma-1)/sigma*psi);

epsw = (1-alpha)/ay_ratio*(1-tauH)/(1+1/(eps_H-(sigma-1)/sigma*psi)); 

kai = psi*(1-sigma)/(sigma*eps_H+psi*(1-sigma));

gamma1 = ((1-kai)*(1-beta_tilda)/epsc-beta_tilda/sigma)*R_tilda;
gamma2 = (1-kai)*(1-beta_tilda)/epsc*(epsw + epsc*kai/(1-kai));


%% endogenous variables

% jump    %%PK: ordering changed

rk = 1; % rental rate 
wage = 2;
output = 3;
gt = 4;  % government expenditure
hours = 5;
cons = 6;
invst = 7;
caput = 8; % capacity utilization

% state
cap = 9;


%% exogenous variables
gamma_x = 10;  % productivity shock

% bonds

bt = 11;   % bond holding

% government expenditure shocks

%gamma_g = 14;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load dynare RE solution here
% Ordering of variables in Dynare file 
% 1   2    3  4    5  6  7   8  9   10   11   12   13  14  15
% y   inv  w  y_H  r  U  a   g  gy  by    k    b   c   H    rk
% Ordering here:
% 1  2    3      4   5    6    7      8     9    10      11   
% rk wage output g  hours cons invst  caput cap  gamma_x b   
%  
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 load coeff_k; 
 coeff_vec_k_RE = zeros(11,2);   % coefficient on capital and bonds
 %coeff_vec_b_RE = zeros(14,2);   % coefficient on bonds
 coeff_vec_s_RE = zeros(11,1);   % coefficient on productivity shocks
coeff_vec_g_RE = zeros(11,1);   % coefficient on government purchase shocks

 
 % 1*2 RE coefficient vector
 coeff_vec_k_RE(1,:) = coeff_k(15,:);   %  rk 1*2 vector
 coeff_vec_k_RE(2,:) = coeff_k(3,:);    %  wage
 coeff_vec_k_RE(3,:) = coeff_k(1,:);    % output
 coeff_vec_k_RE(4,:) = coeff_k(8,:);    %  government expenditure
 coeff_vec_k_RE(5,:) = coeff_k(14,:);   % hours
 coeff_vec_k_RE(6,:) = coeff_k(13,:);    % cons 
 coeff_vec_k_RE(7,:) = coeff_k(2,:);    %  invst 
 coeff_vec_k_RE(8,:) = coeff_k(6,:);    % caput
 coeff_vec_k_RE(9,:) = coeff_k(11,:);  % cap
 coeff_vec_k_RE(10,:) = [0 0] ;    %  gamma_x
 coeff_vec_k_RE(11,:) = coeff_k(12,:);    % g  
  
 % I.i.d government expenditure shocks are in the RE model but not the
 % learning model
 
 coeff_vec_s_RE(1) = coeff_u(15,1);   %  rk 
 coeff_vec_s_RE(2) = coeff_u(3,1);    %  wage
 coeff_vec_s_RE(3) = coeff_u(1,1);    % output
 coeff_vec_s_RE(4) = coeff_u(8,1);    %  gt
 coeff_vec_s_RE(5) = coeff_u(14,1);   % hours
 coeff_vec_s_RE(6) = coeff_u(13,1);    % cons 
 coeff_vec_s_RE(7) = coeff_u(2,1);    %  invst 
 coeff_vec_s_RE(8) = coeff_u(6,1);    % caput
 coeff_vec_s_RE(9) = coeff_u(11,1);  % cap
 coeff_vec_s_RE(10) = 1;    %  gamma_x
 coeff_vec_s_RE(11) = coeff_u(12,1);    % b 
 
  
 coeff_vec_g_RE(1) = coeff_u(15,2);   %  rk 
 coeff_vec_g_RE(2) = coeff_u(3,2);    %  wage
 coeff_vec_g_RE(3) = coeff_u(1,2);    % output
 coeff_vec_g_RE(4) = coeff_u(8,2);    %  gt
 coeff_vec_g_RE(5) = coeff_u(14,2);   % hours
 coeff_vec_g_RE(6) = coeff_u(13,2);    % cons 
 coeff_vec_g_RE(7) = coeff_u(2,2);    %  invst 
 coeff_vec_g_RE(8) = coeff_u(6,2);    % caput
 coeff_vec_g_RE(9) = coeff_u(11,2);  % cap
 %coeff_vec_g_RE(10) = 1;    %  gamma_g
 coeff_vec_g_RE(11) = coeff_u(12,2);    % b 

 %coeff_vec_s_RE(14,:) = [0 1];    % 

 
% Calculating asymptotic moment matrix

% phi matrix 
ini_wR = [(1-coeff_vec_k_RE(2,1)-coeff_vec_k_RE(2,2))*log(gamma); ...
    (-coeff_vec_k_RE(1,1)-coeff_vec_k_RE(1,2))*log(gamma)]...
    + [coeff_vec_k_RE(2,1) coeff_vec_k_RE(2,2); coeff_vec_k_RE(1,1) coeff_vec_k_RE(1,2)]...
    *inv([1-coeff_vec_k_RE(9,1)  -coeff_vec_k_RE(9,2); -coeff_vec_k_RE(11,1) 1-coeff_vec_k_RE(11,2)])...
    *[(1-coeff_vec_k_RE(9,1)-coeff_vec_k_RE(9,2))*log(gamma); (1-coeff_vec_k_RE(11,1)-coeff_vec_k_RE(11,2))*log(gamma)];


phimatrix_RE = [coeff_vec_k_RE(9,:); coeff_vec_k_RE(11,:)];  % see pp. 4 of 22/02/2017 note
omegamatrix_RE = [(1+coeff_vec_s_RE(9))^2 (1+coeff_vec_s_RE(11))*(1+coeff_vec_s_RE(9));...
    (1+coeff_vec_s_RE(11))*(1+coeff_vec_s_RE(9)) (1+coeff_vec_s_RE(11))^2]*sigma_x^2;
% see pp. 4 of 22/02/2017 note
cvec_RE = [(1-coeff_vec_k_RE(9,1)-coeff_vec_k_RE(9,2))*log(gamma);-coeff_vec_k_RE(11,1)*log(gamma)+(1-coeff_vec_k_RE(11,2))*log(gamma)];

%muvec_RE = inv(eye(2)-phimatrix_RE)*cvec_RE; % vector of the mean of capital and bond growth
muvec_RE = [log(gamma);log(gamma)];
vecsigma_RE = inv(eye(4)-kron(phimatrix_RE, phimatrix_RE))* reshape(omegamatrix_RE, [4,1]); 

% asymptotic moment matrix formula
mom_mat_RE = [1 muvec_RE(1) muvec_RE(2); ...
    muvec_RE(1) vecsigma_RE(1)+muvec_RE(1)^2 vecsigma_RE(2)+muvec_RE(1)*muvec_RE(2) ; ...
    muvec_RE(2) vecsigma_RE(2)+muvec_RE(1)*muvec_RE(2) vecsigma_RE(4)+muvec_RE(2)^2];

%  RE coefficient in the law of motion for capital and bond growth
coeff_mat_k_AL_ini =[cvec_RE(1); phimatrix_RE(1,:)']; 

coeff_mat_b_AL_ini =[cvec_RE(2);phimatrix_RE(2,:)']; 

%  RE coefficient in the law of motion for wage and rental rate growth

coeff_mat_w_AL_ini =[(1-coeff_vec_k_RE(2,1)-coeff_vec_k_RE(2,2))*log(gamma);coeff_vec_k_RE(2,1);coeff_vec_k_RE(2,2)];

coeff_mat_R_AL_ini =[(-coeff_vec_k_RE(1,1)-coeff_vec_k_RE(1,2))*log(gamma);coeff_vec_k_RE(1,1);coeff_vec_k_RE(1,2)];  %rk

% output & productivity
coeff_mat_y_AL_ini = [(1-coeff_vec_k_RE(3,1)-coeff_vec_k_RE(3,2))*log(gamma);coeff_vec_k_RE(3,1);coeff_vec_k_RE(3,2)];
coeff_mat_pr_AL_ini = [(1-coeff_vec_k_RE(2,1)-coeff_vec_k_RE(2,2))*log(gamma);coeff_vec_k_RE(2,1);coeff_vec_k_RE(2,2)];

% initial moments: asymptotic moment matrix
mom_control = 2;
if mom_control == 1
  Mom_mat_ini = eye(3);     
elseif  mom_control == 2 
   Mom_mat_ini  = mom_mat_RE; 
end



%% GENERATES SHOCKS
%Note data are simulated for long periods to calculate the output
%volatility

randn('state',9)

n1 = 500; % no. of repetitions

 nsim =2162;

 epsZ = randn(n1,nsim); 
 shockmatx = zeros(n1,nsim); 
  shockmat = epsZ*sigma_x; % the shock sequences
  
Tend = nsim; 
%lastT = 162; % using the lastT periods of data to calculate the statistics
%lastT = 162;
lastT = 162;
Tend1 = Tend-lastT+1;    % The starting quarter of shorted data for calculating statistics

% gain parameters 

gain_control =1  ;% 1: constant gain 0 : decreasing gain
if gain_control ==1 
    
% 
%  gaink = 0.014; 
%  gainb = 0.014; 
%  gainw = 0.014;% 0.011;%0.004; %0
%  gainR = 0.04;% 0.011;%0.005; %0
%  gainkm = 0.014; 
%  gainy = 0.014; %  does not matter for the temporary equilibrium
%   gainy1 = gainw +0.019; 
%  gaink = 0.008; 
%  gainb = 0.008; 
%  gainw = 0.008;% 0.011;%0.004; %0
%  gainR = 0.03;% 0.011;%0.005; %0
%  gainkm = 0.008; 
%  gainy = 0.008; %  does not matter for the temporary equilibrium
%   gainy1 = gainw +0.028; 


 gaink = 0.003; 
 gainb = 0.003; 
 gainw = 0.003;% 0.011;%0.004; %0
 gainR = 0.015;% 0.011;%0.005; %0
 gainkm = 0.003; 
 gainy = 0.003; %  does not matter for the temporary equilibrium
 gainy1 = gainw +0.032; 


else 
   gaink = 1./(200+(1:nsim)); 
   gainb = 1./(200+(1:nsim)); 
   gainw =  1./(200+(1:nsim));
       gainR =  1./(200+(1:nsim));
          gainkm =  1./(200+(1:nsim));
    
end

 %Generating IRFs
 irfcontrol = 1;        % 1: generating IRFs for three models: RE, 
                         %Learning and only underparameterization
 shockcontrol = 7; 
 if irfcontrol==1 
     

%  xx = [0 0 0 0 0];
%  main_calibration_km_IRF(xx);
 
 xx = [gainw gainR gaink gainkm gainb];  
  if shockcontrol == 1   %   IRF following one-off shock
      main_calibration_km_IRF(xx);
  elseif shockcontrol == 2 % IRF following a sequence of shock 
      main_calibration_km_IRF_great_recession(xx);
  elseif  shockcontrol == 3                  % IRF following a sequence of shocks and fiscal consolidation
       main_calibration_km_IRF_great_recession_fc(xx);
  elseif shockcontrol ==4
      main_calibration_km_IRF_great_recession_heterogeneous(xx);
  elseif shockcontrol ==5;
      main_calibration_km_IRF_change_g(xx);
  elseif shockcontrol ==6;
      main_calibration_km_IRF_severity(xx);      
  elseif shockcontrol ==7;
      main_calibration_km_IRF_excessive(xx);
  elseif shockcontrol ==8;
      main_calibration_km_IRF_excessive2(xx); 
  elseif shockcontrol ==9;
      
%      xx = [xx gainy1]; 
      main_calibration_km_optimal_gain(xx); 
     
      
  end
 else 
  
   
 end
