clc;
clear all;
close all;

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

load coeff_k;
load coeff_u; 
load coeff_mat_k_AL2000.mat;


load coeff_mat_b_AL2000.mat;
load coeff_mat_R_AL2000.mat;
load coeff_mat_w_AL2000.mat;

load y_tg2000.mat; 
load Dlogy2000.mat; 
load Dlogw2000.mat;
load DlogR2000.mat;
load Dlogk1999.mat;
load Dlogb1999.mat;
load Dlogk2000.mat;
load Dlogb2000.mat;

load  khat2000.mat ;
load  bhat2000.mat ;
load what2000.mat;
load RKhat2000.mat;
load yhat2000.mat;



delta = 0.025; %% depreciation rate

alpha = 0.38;  %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.4;

%tauH = 0.4;
tauH = 0.42;

tauC= 0.17;  % 0.17
gy_ratio = 0.18; 
 
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.75;%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

 
 % 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_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 = 1;
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',10)

n1 = 500; % no. of repetitions

 nsim =540;
NN= 40; % evaluation span
 epsZ = randn(n1,nsim); 
 shockmatx = zeros(n1,nsim); 
  shockmat = epsZ*sigma_x; % the shock sequences
  std_A = sigma_x;

  
ngain =100; 
og_bias1 = zeros(n1,nsim);  % bias in output gap for an individual value of the gain parameter
og_bias_ave = zeros(ngain, 1); % mean bias
psiy_ec = zeros(n1,nsim);  %  
gain_grids= linspace(0, 0.05, ngain);

% gain parameters 

%  gaink = 0.007; 
%  gainb = 0.007; 
%  gainw = 0.007;% 0.011;%0.004; %0
%  gainR = 0.016;% 0.011;%0.005; %0
%  gainkm = 0.007; 
%  gainy = 0.007; %  does not matter for the temporary equilibrium
 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; 
 
 for j = 1: ngain 
   gainy1 = gain_grids(j); %  does not matter for the temporary equilibrium
%   gaink = 0.005; 
%  gainb = 0.005; 
%  gainw = 0.005;% 0.011;%0.004; %0
%  gainR = 0.01;% 0.011;%0.005; %0
%  gainkm = 0.005; 
%  gainy1 = gainw + 0.03; %  does not matter for the temporary equilibrium
% gaink = 0.008; 
%  gainb = 0.008; 
%  gainw = 0.008;% 0.011;%0.004; %0
%  gainR = 0.0016;% 0.011;%0.005; %0
%  gainkm = 0.008; 
%  gainy = 0.008; %  does not matter for the temporary equilibrium



%

%   khat1 = 0;  % intial capital stock
%   ahat1 = 0; % initial asset holdings 
%   bhat1 = 0 ; %initial bond stocks

% note we need two values for wage and interst rates to initialize the
% learning algorithm
% 
% wmin1 = 0;  % wage at period -1
% w0 = 0;  % wage at period 0 
% 
% rmin1 = 0 ;  % interest rate at period -1
% r0 = 0;  % interest rate at period 0
% 
% kmin1=0;
% k0 =0;
% 
% shock0 = 0; % intial shock value 

%%Beliefs : Learning model
% psiw = zeros(n1,nsim);  % record agents' wage growth belief
% psir = zeros(n1,nsim);  % record agents' interest rate growth belief
coeff_mat_w_AL = zeros(3, nsim);
coeff_mat_R_AL = zeros(3, nsim);
coeff_mat_k_AL = zeros(3, nsim);
coeff_mat_b_AL = zeros(3, nsim);


%Mom_mat_k = zeros(2,2*nsim);

% Record DlogR, DlogW, the dependent variables

Dlogw = zeros(n1, nsim);   



DlogR = zeros(n1, nsim);

Dlogk = zeros(n1, nsim);



Dlogb = zeros(n1, nsim);


%regressors in the learning model (income equation): set to the steady
%state

 
%%Endogenous variables : learning model

% record simulated series: number of repetitions * no. of simul periods
%
% consumption, output, interst rate, hours, wage, rental rate, investment,
% infinite sums of y, r, capital stock, productivity




y_tg = zeros(n1,nsim);
xtg = zeros(n1, nsim);


RKhat = zeros(n1,nsim);
what = zeros(n1,nsim);
yhat = zeros(n1,nsim);
ghat = zeros(n1,nsim);
Hhat = zeros(n1,nsim);
chat = zeros(n1,nsim);    
ihat = zeros(n1,nsim);
caputhat = zeros(n1,nsim);
khat = zeros(n1,nsim);
Prhat = zeros(n1,nsim);
Swe = zeros(n1,nsim);
Sre = zeros(n1,nsim);
bhat = zeros(n1, nsim); 


for i = 1:n1 % loop over repetions
    Dlogw0 = Dlogw2000(i);
    DlogR0 = DlogR2000(i); 
    Dlogk0 = Dlogk2000(i);
    Dlogkm1 = Dlogk1999(i); 

    Dlogb0 = Dlogb2000(i);
    Dlogbm1 = Dlogb1999(i); 
    y_tg0 = y_tg2000(i);
    x_vec_m1 = [1;Dlogkm1; Dlogbm1]; 
    x_vec_ini = [1; Dlogk0; Dlogb0];

    for t = 1:nsim  % loop over simul periods
       if t==1 
  %         khat(i,t) = khat1; 
           % belief updating in wage equation
 %          Mom_mat_k(:, 1 + 2*(t-1):2*t) = Mom_mat_k_ini + gainkm *(x_vec_m1*x_vec_m1'-Mom_mat_k_ini);
          
           coeff_mat_k_AL(:,t) = coeff_mat_k_AL2000(i,:)' + gaink*x_vec_m1*(Dlogk0-coeff_mat_k_AL2000(i,:)*x_vec_m1);
           coeff_mat_b_AL(:,t) = coeff_mat_b_AL2000(i,:)' + gainb*x_vec_m1*(Dlogb0-coeff_mat_b_AL2000(i,:)*x_vec_m1);
           % belief updating in interest rate equation
           coeff_mat_R_AL(:,t) = coeff_mat_R_AL2000(i,:)' + gainR*x_vec_m1*(DlogR0-coeff_mat_R_AL2000(i,:)*x_vec_m1) ;
           
           coeff_mat_w_AL(:,t) = coeff_mat_w_AL2000(i,:)' + gainw*x_vec_m1*(Dlogw0-coeff_mat_w_AL2000(i,:)*x_vec_m1) ;
           

          y_tg(i,t) = y_tg0+gainy1 * (Dlogy2000(i)-y_tg0); 

       elseif t==2
          
   %        Mom_mat_k(:, 1 + 2*(t-1):2*t) = Mom_mat_k(:, 1 + 2*(t-2):2*(t-1)) + gainkm*(x_vec_ini*x_vec_ini' - Mom_mat_k(:, 1 + 2*(t-2):2*(t-1)));
           
           coeff_mat_k_AL(:,t) = coeff_mat_k_AL(:,t-1) + gaink*x_vec_ini*(Dlogk(i,t-1)-coeff_mat_k_AL(:,t-1)'*x_vec_ini);
           coeff_mat_b_AL(:,t) = coeff_mat_b_AL(:,t-1) + gainb*x_vec_ini*(Dlogb(i,t-1)-coeff_mat_b_AL(:,t-1)'*x_vec_ini);
           % belief updating in interest rate equation
           coeff_mat_R_AL(:,t) = coeff_mat_R_AL(:,t-1)+ gainR*x_vec_ini*(DlogR(i,t-1)-coeff_mat_R_AL(:,t-1)'*x_vec_ini) ;
           
           coeff_mat_w_AL(:,t) = coeff_mat_w_AL(:,t-1) + gainw*x_vec_ini*(Dlogw(i,t-1)-coeff_mat_w_AL(:,t-1)'*x_vec_ini);
           
           y_tg(i,t) = y_tg(i,t-1)+gainy1*(Dlogy(i,t-1)-y_tg(i,t-1)); 

       else
%            if coeff_mat_k_AL(:,t-1) == coeff_mat_k_AL(:,t-2)
%                x_vec = [1; Dlogk(i,t-3)]; 
%            else
           x_vec = [1; Dlogk(i,t-2); Dlogb(i,t-2)]; 
%            end
           
     %      Mom_mat_k(:, 1 + 2*(t-1):2*t) = Mom_mat_k(:, 1 + 2*(t-2):2*(t-1)) + gainkm*(x_vec*x_vec' - Mom_mat_k(:, 1 + 2*(t-2):2*(t-1)));
%            Mom_mat_k(:, 1 + 2*(t-1):2*t)
           coeff_mat_k_AL(:,t) = coeff_mat_k_AL(:,t-1) + gaink*x_vec*(Dlogk(i,t-1)-coeff_mat_k_AL(:,t-1)'*x_vec);
           coeff_mat_b_AL(:,t) = coeff_mat_b_AL(:,t-1) + gainb*x_vec*(Dlogb(i,t-1)-coeff_mat_b_AL(:,t-1)'*x_vec);
                       
           % belief updating in interest rate equation
           coeff_mat_R_AL(:,t) = coeff_mat_R_AL(:,t-1)+ gainR*x_vec*(DlogR(i,t-1)-coeff_mat_R_AL(:,t-1)'*x_vec);
           
           coeff_mat_w_AL(:,t) = coeff_mat_w_AL(:,t-1) + gainw*x_vec*(Dlogw(i,t-1)-coeff_mat_w_AL(:,t-1)'*x_vec);
            y_tg(i,t) = y_tg(i,t-1)+gainy1*(Dlogy(i,t-1)-y_tg(i,t-1)); 

       end
       
       %%%Trend growth forecasts
       

           xw = coeff_mat_w_AL(:,t);  % wage growth belief
           xR = coeff_mat_R_AL(:,t);  % interest growth belief
           xk = coeff_mat_k_AL(:,t); 
           xb = coeff_mat_b_AL(:,t); 
           if t<=500 
             xtg(i,t) = (t+16)*(y_tg(i,t)-log(gamma)); 
           elseif t>500
             xtg(i,t) = (t-500+16)*(y_tg(i,t)-log(gamma));   
           end
           if t==1 
               ini_con = [shockmat(i,t) khat2000(i) bhat2000(i) xtg(i,t)]; 
           else
               ini_con = [shockmat(i,t) khat(i,t-1) bhat(i,t-1) xtg(i,t)]; 
           end

%            end
           %      xw : trend belief in w equation
           %      xR : trend belief in R equation
           % ini_con: ordering gam lagk
           [y_vec1, fore_dis] = Learningsol_AL(xw, xR, xk, xb, ini_con); 
           % solutions
           %Note: ordering of the variables: chat, yhat, rhat, Hhat, what, RKhat, ihat
           %y_vec1
           RKhat(i,t) = y_vec1(1);
           what(i,t) = y_vec1(2);
           yhat(i,t) = y_vec1(3);
           ghat(i,t) = y_vec1(4);
           Hhat(i,t) = y_vec1(5);
           chat(i,t) = y_vec1(6);  
           ihat(i,t) = y_vec1(7);
           caputhat(i,t) = y_vec1(8);
           khat(i,t) = y_vec1(9);
           bhat(i,t) = y_vec1(11);
           
           swe(i,t) = fore_dis(1);
           sre(i,t) = fore_dis(2) ; 
           
          if t==1
              
           Dlogw(i,t) = what(i,t)-what2000(i)+log(gamma)+shockmat(i,t); 
           DlogR(i,t) = RKhat(i,t)-RKhat2000(i); 
           Dlogk(i,t) = khat(i,t)-khat2000(i)+log(gamma)+shockmat(i,t);   
           Dlogy(i,t) = yhat(i,t) - yhat2000(i) + log(gamma)+shockmat(i,t);   
           Dlogb(i,t) = bhat(i,t) - bhat2000(i) + log(gamma)+shockmat(i,t);   
          elseif   t>=2
           Dlogw(i,t) = what(i,t)-what(i,t-1)+log(gamma)+shockmat(i,t); 
           DlogR(i,t) = RKhat(i,t)-RKhat(i,t-1); 
           Dlogk(i,t) = khat(i,t)-khat(i,t-1)+log(gamma)+shockmat(i,t); 
           Dlogy(i,t) = yhat(i,t) - yhat(i,t-1) + log(gamma)+shockmat(i,t);           
           Dlogb(i,t) = bhat(i,t) - bhat(i,t-1) + log(gamma)+shockmat(i,t); 
          end


    end % end Tend loop over simul periods

    
    
end % end n1 loop over repetition

         og_bias_ave(j) =sum(reshape(abs(xtg(:,end-NN+1:end)), [1, n1*NN]))/n1/NN;

 end

figure 
set(gca,'fontsize',18)
h0= plot(gain_grids, og_bias_ave)
set(h0,'linewidth',1.8)
%title('Mean absolute bias')
h1=xlabel('Gain parameter')
set(h1,'fontsize',18)
h2= ylabel('Mean absolute estimation error')
set(h2,'fontsize',18)
[m,index] = min(og_bias_ave); 
%saveas(gcf, 'optimal_gain.eps','epsc')

display('optimal gain minimizing the mean absolute bias')
gain_grids(index)
display('Mean absolute bias associated with the optimal gain')

saveas(gcf,'optimalgain','epsc')
