
%created 21 February 2017: 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function main_calibration_km_IRF(xx)

global coeff_vec_k_RE coeff_vec_s_RE muvec_RE  gamma ...
 coeff_mat_k_AL_ini coeff_mat_b_AL_ini   coeff_mat_R_AL_ini   coeff_mat_w_AL_ini mom_mat_RE

gainw = xx(1);
gainR = xx(2);
gaink = xx(3);
gainkm = xx(4); 
gainb = xx(5) ; 
gainy1 = gainw+0.03;

Font_size = 11;  
Line_width = 1.5; 

%% GENERATES SHOCKS
n1 =1;
%randn('state',9)
 nsim =40;
 %epsZ = randn(n1,nsim); 
  Tend = nsim;

shockmat = zeros(n1,nsim); 

shockmat(1) = -0.003% std_A;
% shockmat(2) = -0.003;
% shockmat(3) = -0.003;
% shockmat(4) = -0.003;
% shockmat(5) = -0.003;
% shockmat(6) = -0.003;
% shockmat(7) = -0.003;
% shockmat(8) = -0.007;


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


% initial beliefs


coeff_mat_k_AL = zeros(3, nsim);
coeff_mat_b_AL = zeros(3, nsim); 


coeff_mat_w_AL = zeros(3, nsim);
coeff_mat_R_AL = zeros(3, nsim);

mom_mat_ini = mom_mat_RE; 
mom_mat_kb = zeros(3, 3*nsim); 

% Record DlogR, DlogW, the dependent variables

Dlogw = zeros(1, nsim); 
Dlogw0 = log(gamma);

DlogR = zeros(1, nsim);
DlogR0 = 0; 


Dlogk = zeros(1, nsim);
Dlogk0 = muvec_RE(1);
Dlogkm1 = muvec_RE(1); % in period minus one

Dlogb = zeros(1,nsim);
Dlogb0 = muvec_RE(2);
Dlogbm1 = muvec_RE(2); 

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

x_vec_m1 = [1;Dlogkm1; Dlogbm1];  % regressors in period minus one 
x_vec_ini = [1; Dlogk0 ; Dlogb0];  % regressors in period zero

y_tg0 = log(gamma);  % initial belief about trend growth
% y_tg = [log(gamma)*ones(1,16) zeros(n1,nsim)];  % belief vector about trend growth
y_tg = zeros(n1,nsim);  % belief vector about trend growth

RKhat = zeros(n1,nsim);  %  rental rate
what = zeros(n1,nsim);   % wages 
yhat = zeros(n1,nsim);   %  output
%rthat = zeros(n1,nsim);   %interest rate
Hhat = zeros(n1,nsim);   % hours 
chat = zeros(n1,nsim);   % consumption    
ihat = zeros(n1,nsim);   % investment
caputhat = zeros(n1,nsim); % capacity utilization
khat = zeros(n1,nsim);     % capital 
ghat = zeros(n1,nsim);     % government spending
%ahat = zeros(n1,nsim);     % assets
bhat = zeros(n1,nsim);     % bond holdings

Prhat = zeros(n1,nsim);    %  productivity
Swe = zeros(n1,nsim);      %  forecast of discounted sum of wage rates
Sre = zeros(n1,nsim);      %  forecast of discounted sum of rental rates

y_trend_o = zeros(n1,nsim+16);
y_trend_s = zeros(n1, nsim+16);
xtg = zeros(n1,nsim);

    
    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);
%          size(x_vec_m1)
%          size(coeff_mat_k_AL_ini)
           coeff_mat_k_AL(:,t) = coeff_mat_k_AL_ini + gaink*x_vec_m1*(Dlogk0-coeff_mat_k_AL_ini'*x_vec_m1);
           coeff_mat_b_AL(:,t) = coeff_mat_b_AL_ini + gainb*x_vec_m1*(Dlogb0-coeff_mat_b_AL_ini'*x_vec_m1);
           
           mom_mat_kb(:, 1+3*(t-1):3*t) = mom_mat_ini;
           % belief updating in interest rate equation
           coeff_mat_R_AL(:,t) = coeff_mat_R_AL_ini + gainR*x_vec_m1*(DlogR0-coeff_mat_R_AL_ini'*x_vec_m1) ;
           
           coeff_mat_w_AL(:,t) = coeff_mat_w_AL_ini + gainw*x_vec_m1*(Dlogw0-coeff_mat_w_AL_ini'*x_vec_m1) ;
           y_tg(t) = y_tg0; 
           

           
       elseif t==2
           mom_mat_kb(:, 1+3*(t-1):3*t) = mom_mat_ini;
           
           %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(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(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(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(t-1)-coeff_mat_w_AL(:,t-1)'*x_vec_ini);
           y_tg(t) = y_tg(t-1) +gainy1*(DlogY(t-1)-y_tg(t-1)); 
           

       else
           x_vec = [1; Dlogk(t-2); Dlogb(t-2)]; 
           mom_mat_kb(:, 1+3*(t-1):3*t) = mom_mat_kb(:, 1+3*(t-2):3*(t-1)) + gainkm * (x_vec * x_vec' - mom_mat_kb(:, 1+3*(t-2):3*(t-1)));
        %   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)));
           
           coeff_mat_k_AL(:,t) = coeff_mat_k_AL(:,t-1) + gaink*x_vec*(Dlogk(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(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(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(t-1)-coeff_mat_w_AL(:,t-1)'*x_vec);
           y_tg(t) = y_tg(t-1) +gainy1*(DlogY(t-1)-y_tg(t-1)); 
       end      
           
       % beliefs about parameters in PLM for the growth of wage, bond holding, rental rates and capital  
           xw = coeff_mat_w_AL(:,t);  % wage growth 
           xb = coeff_mat_b_AL(:,t) ; % bond holding growth 
           xR = coeff_mat_R_AL(:,t)  ;% interest growth 
           xk = coeff_mat_k_AL(:,t)  ;% capital growth
           xtg(t) = (t+16)*(y_tg(t)-log(gamma)); % feed in 4 years before crisis
           if t==1
              ini_con = [shockmat(t) khat1 bhat1 xtg(t)] ; 
           else
              ini_con = [shockmat(t) khat(t-1) bhat(t-1) xtg(t)]; 
%                display('learning model ahat(i,t-1) bhat(i,t-1)')
%            ahat(i,t-1) 
%            bhat(i,t-1)
           end
           %      xw : trend belief in w equation
           %      xR : trend belief in R equation
           % ini_con: ordering gam lagk
%            display('ini_con')
%            format long
%            ini_con
           [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
%            display('forecast of wage sums')
%            fore_dis(1)
%            display('forecast of rental sums')
%            fore_dis(2)
%            
%        display('y_vec1')
%        y_vec1
%            
           
           RKhat(t) = y_vec1(1);
           what(t) = y_vec1(2);
           yhat(t) = y_vec1(3);
           ghat(t) = y_vec1(4);
           Hhat(t) = y_vec1(5);
           chat(t) = y_vec1(6);  
           ihat(t) = y_vec1(7);
           caputhat(t) = y_vec1(8);
           khat(t) = y_vec1(9);
           bhat(t) = y_vec1(11);

           
           Prhat(t) = yhat(t)-Hhat(t); 
           Swe(t) = fore_dis(1); 
           Sre(t) = fore_dis(2); 
           
           
    
           
          if t==1
           DlogY(t) = yhat(t)-0+log(gamma)+shockmat(t); 
           Dlogw(t) = what(t)-0+log(gamma)+shockmat(t); 
           DlogR(t) = RKhat(t)-0; 
           Dlogk(t) = khat(t)-0+log(gamma)+shockmat(t); 
           Dlogb(t) = bhat(t)-0+log(gamma)+shockmat(t); 
          elseif   t>=2
           DlogY(t) = yhat(t)-yhat(t-1)+log(gamma)+shockmat(t); 
           Dlogw(t) = what(t)-what(t-1)+log(gamma)+shockmat(t); 
           DlogR(t) = RKhat(t)-RKhat(t-1); 
           Dlogk(t) = khat(t)-khat(t-1)+log(gamma)+shockmat(t);
           Dlogb(t) = bhat(t)-bhat(t-1)+log(gamma)+shockmat(t);
          end

    

    end % end Tend loop over simul periods

     [y_tg_min , index] = min(y_tg);
%      display('index')
%      index
    for k = 1: t+16
       y_trend_s(k) = k*y_tg(index);
       y_trend_o(k) = k*log(gamma); 
    end
    for k =1:16
        xtg_past(k) = k*(y_tg(index)-log(gamma)); 
    end

% display('coeff_mat_k_AL')
% coeff_mat_k_AL
% 
% display('coeff_mat_b_AL')
% coeff_mat_b_AL
% 
% display('coeff_mat_R_AL')
% coeff_mat_R_AL
% 
% display('coeff_mat_w_AL')
% coeff_mat_w_AL
% %Initial conditions
% 
display('y_tg')
y_tg

display('RE model@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@')
HRE = zeros(n1, Tend);
yyRE  = zeros(n1, Tend);
cRE  = zeros(n1, Tend);
wRE = zeros(n1, Tend);
REinv = zeros(n1, Tend);
PrRE = zeros(n1, Tend);
kRE = zeros(n1,Tend);
caputRE = zeros(n1,Tend);
RkRE = zeros(n1,Tend);
%RtRE = zeros(n1,Tend);
gtRE = zeros(n1,Tend);
%atRE = zeros(n1,Tend);
btRE = zeros(n1,Tend);

SwRE = zeros(n1,Tend);
SrRE = zeros(n1,Tend);




RKhat_RE2 = zeros(n1,nsim);  %  rental rate
what_RE2 = zeros(n1,nsim);   % wages 
yhat_RE2 = zeros(n1,nsim);   %  output
%rthat = zeros(n1,nsim);   %interest rate
Hhat_RE2 = zeros(n1,nsim);   % hours 
chat_RE2 = zeros(n1,nsim);   % consumption    
ihat_RE2 = zeros(n1,nsim);   % investment
caputhat_RE2 = zeros(n1,nsim); % capacity utilization
khat_RE2 = zeros(n1,nsim);     % capital 
ghat_RE2 = zeros(n1,nsim);     % government spending
%ahat = zeros(n1,nsim);     % assets
bhat_RE2 = zeros(n1,nsim);     % bond holdings



% xw_RE = [(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)];
% xR_RE = [(-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)]; 
% xk_RE =  [cvec_RE(1); phimatrix_RE(1,:)'];
% xb_RE = [cvec_RE(2);phimatrix_RE(2,:)']; 

xw_RE = coeff_mat_w_AL_ini
xR_RE = coeff_mat_R_AL_ini;
xk_RE = coeff_mat_k_AL_ini;
xb_RE = coeff_mat_b_AL_ini;

for j = 1:nsim
   
       if j==1    
            HRE(j) = coeff_vec_s_RE(5)*shockmat(j);
            yyRE(j) = coeff_vec_s_RE(3)*shockmat(j) ; 
            cRE(j) =coeff_vec_s_RE(6)*shockmat(j);
            wRE(j) =coeff_vec_s_RE(2)*shockmat(j);
            caputRE(j) = coeff_vec_s_RE(8)*shockmat(j); 
      %      RtRE(j,i) = coeff_vec_s_RE(4)*shockmat(j,i);
            REinv(j) = coeff_vec_s_RE(7)*shockmat(j);
            kRE(j) = coeff_vec_s_RE(9)*shockmat(j);
            gtRE(j) = coeff_vec_s_RE(4)*shockmat(j);
      %      atRE(j,i) = coeff_vec_s_RE(12)*shockmat(j,i);
            btRE(j) = coeff_vec_s_RE(11)*shockmat(j);
            RkRE(j) = coeff_vec_s_RE(1)*shockmat(j);

                ini_con_RE = [shockmat(j) 0  0]; 

               [y_vec1_RE, fore_dis_RE]  = Learningsol(xw_RE, xR_RE, xk_RE, xb_RE, ini_con_RE);

      
          
          
               SwRE(j) = fore_dis_RE(1);
               
               SrRE(j) = fore_dis_RE(2);                
            
       else
            HRE(j) = coeff_vec_k_RE(5,1)*kRE(j-1) + coeff_vec_k_RE(5,2)*btRE(j-1)+ coeff_vec_s_RE(5)*shockmat(j);
            yyRE(j) = coeff_vec_k_RE(3,1)*kRE(j-1) + coeff_vec_k_RE(3,2)*btRE(j-1) + coeff_vec_s_RE(3)*shockmat(j);          
            cRE(j) = coeff_vec_k_RE(6,1)*kRE(j-1) + coeff_vec_k_RE(6,2)*btRE(j-1) + coeff_vec_s_RE(6)*shockmat(j);         
            wRE(j) = coeff_vec_k_RE(2,1)*kRE(j-1)+coeff_vec_k_RE(2,2)*btRE(j-1) +coeff_vec_s_RE(2)*shockmat(j);                      
            REinv(j) = coeff_vec_k_RE(7,1)*kRE(j-1)+coeff_vec_k_RE(7,2)*btRE(j-1) + coeff_vec_s_RE(7)*shockmat(j);                      
            kRE(j) = coeff_vec_k_RE(9,1)*kRE(j-1)+ coeff_vec_k_RE(9,2)*btRE(j-1) +coeff_vec_s_RE(9)*shockmat(j); 
            gtRE(j) = coeff_vec_k_RE(4,1)*kRE(j-1) + coeff_vec_k_RE(4,2)*btRE(j-1)+ coeff_vec_s_RE(4)*shockmat(j);
         %   atRE(j,i) =  coeff_vec_k_RE(12,1)*kRE(j,i-1) + coeff_vec_k_RE(12,2)*btRE(j,i-1)+ coeff_vec_s_RE(12)*shockmat(j,i);
            btRE(j) =  coeff_vec_k_RE(11,1)*kRE(j-1) + coeff_vec_k_RE(11,2)*btRE(j-1)+ coeff_vec_s_RE(11)*shockmat(j);
            RkRE(j) = coeff_vec_k_RE(1,1)*kRE(j-1)+coeff_vec_k_RE(1,2)*btRE(j-1) + coeff_vec_s_RE(1)*shockmat(j);
   
            
               ini_con_RE = [shockmat(j) kRE(j-1) btRE(j-1)]; 
%                display('ini_con_RE')
%                ini_con_RE
               [y_vec1_RE, fore_dis_RE]  = Learningsol(xw_RE, xR_RE, xk_RE, xb_RE, ini_con_RE);
          %@@@@@@@@@@@@ Continue from here @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
          
               SwRE(j) = fore_dis_RE(1);
               
               SrRE(j) = fore_dis_RE(2);         
        
       end
%                  display('shockmat(j)')
%                 shockmat(j)
%        display('y_vec1_RE')
%        y_vec1_RE
           RKhat_RE2(j) = y_vec1_RE(1);
           what_RE2(j) = y_vec1_RE(2);
           yhat_RE2(j) = y_vec1_RE(3);
           ghat_RE2(j) = y_vec1_RE(4);
           Hhat_RE2(j) = y_vec1_RE(5);
           chat_RE2(j) = y_vec1_RE(6);  
           ihat_RE2(j) = y_vec1_RE(7);
           caputhat_RE2(j) = y_vec1_RE(8);
           khat_RE2(j) = y_vec1_RE(9);
           bhat_RE2(j) = y_vec1_RE(11);
   end
   


PrRE =  yyRE  - HRE; 

figure 
subplot(3,3,1)
% capital stock plot

plot(kRE, 'b-o')
hold on
plot(zeros(Tend,1))
hold on

plot(khat,'r-d')
title('Capital')

%legend('RE', 'Learning');
%plot(zeros(Tend,1))

%output plot
subplot(3,3,2)
plot(yyRE, 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(yhat,'r-d')
title('Output')

% investment plot
subplot(3,3,3)
plot(REinv, 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(ihat,'r-d')
title('Investment')

% consumption plot
subplot(3,3,4)
plot(cRE, 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(chat,'r-d')
title('Consumption')
display('++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
display('chat associated with IRF')
display('++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
%format long;
chat(2) 
% hours plot
subplot(3,3,5)
plot(HRE, 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(Hhat,'r-d')
title('Hours')

% wage plot
subplot(3,3,6)
plot(wRE, 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(what,'r-d')
title('Wage')

% rental rate plot
subplot(3,3,7)
plot(RkRE, 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(RKhat,'r-d')
title('Rental Rate')

% income sum plot
subplot(3,3,8)
plot(SwRE, 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(Swe,'r-d')
title('Infinite Sum of wage')

% interest rate sum plot
subplot(3,3,9)
plot(SrRE, 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(Sre,'r-d')
title('Infinite Sum of Interest Rate')

figure 
subplot(3,2,1)
set(gca,'fontsize',Font_size)
h1 = plot(yyRE*100, 'b-o')
set(h1,'linewidth',Line_width)

hold on 
plot(yhat*100,'r-d')
hold on 
plot(zeros(Tend,1))
legend({'RE', 'Learning'}, 'Location', 'Southeast', 'FontSize', 14);
ylhand = get(gca,'ylabel')
set(ylhand,'string','% dev. from SS','fontsize',17) 
title('Output')


% capital stock plot

%plot(zeros(Tend,1))

%output plot

% investment plot
subplot(3,2,2)
set(gca,'fontsize',Font_size)
h2 = plot(REinv*100, 'b-o')
set(h2,'linewidth',Line_width)
hold on 
plot(zeros(Tend,1))
hold on 
plot(ihat*100,'r-d')
title('Investment')

% consumption plot
subplot(3,2,3)
set(gca,'fontsize',Font_size)
h3 = plot(cRE*100, 'b-o')
set(h3,'linewidth',Line_width)
hold on 
plot(zeros(Tend,1))
hold on 
plot(chat*100,'r-d')
title('Consumption')
display('++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
display('chat associated with IRF')
display('++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
%format long;
chat(2) 
% hours plot
subplot(3,2,4)
set(gca,'fontsize',Font_size)
h4 = plot(HRE*100, 'b-o')
set(h4,'linewidth',Line_width)
hold on 
plot(zeros(Tend,1))
hold on 
plot(Hhat*100,'r-d')
title('Hours')




% government expenditures
subplot(3,2,5)
set(gca,'fontsize',Font_size)
hi1 = plot(gtRE*100, 'b-o')
set(hi1,'linewidth',Line_width)

hold on 
plot(ghat*100,'r-d')
hold on 
plot(zeros(Tend,1))
%legend({'RE', 'Learning'}, 'Location', 'Southeast', 'FontSize', 15);

title('Government expenditures')

%plot(zeros(Tend,1))

%output plot

% investment plot
subplot(3,2,6)
set(gca,'fontsize',Font_size)
h2 = plot((btRE-yyRE)*100, 'b-o')
set(h2,'linewidth',Line_width)
hold on 
plot(zeros(Tend,1))
hold on 
plot((bhat-yhat)*100,'r-d')
title('Debt/GDP ratio')
saveas(gcf,'IRF','epsc')
% % consumption plot
% subplot(2,2,3)
% set(gca,'fontsize',Font_size)
% h3 = plot(cRE*100, 'b-o')
% set(h3,'linewidth',Line_width)
% hold on 
% plot(zeros(Tend,1))
% hold on 
% plot(chat*100,'r-d')
% title('Consumption')
% display('++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
% display('chat associated with IRF')
% display('++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
% %format long;
% chat(2) 
% % hours plot
% subplot(2,2,4)
% set(gca,'fontsize',Font_size)
% h4 = plot(HRE*100, 'b-o')
% set(h4,'linewidth',Line_width)
% hold on 
% plot(zeros(Tend,1))
% hold on 
% plot(Hhat*100,'r-d')
% title('Hours')
% saveas(gcf,'IRF','epsc')



figure
% income sum plot
subplot(2,2,1)
plot(SwRE, 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(Swe,'r-d')
legend('RE', 'Learning');
title('Infinite Sum of wage')

% interest rate sum plot
subplot(2,2,2)
plot(SrRE, 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(Sre,'r-d')
title('Infinite Sum of Interest Rate')

figure 
subplot(3,3,1)
% capital stock plot

plot(kRE, 'b-o')
hold on
plot(zeros(Tend,1))
hold on
plot(khat_RE2,'r-d')

display('khat_RE2')
khat_RE2

title('Capital')
%legend('RE', 'Learning');
%plot(zeros(Tend,1))

%output plot
subplot(3,3,2)
plot(yyRE, 'b-o')
hold on 
plot(zeros(Tend,1))
hold on
plot(yhat_RE2,'r-d')

title('Output')

% investment plot
subplot(3,3,3)
plot(REinv, 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(ihat_RE2,'r-d')
title('Investment')

% consumption plot
subplot(3,3,4)
plot(cRE, 'b-o')
hold on 
plot(zeros(Tend,1))
hold on
plot(chat_RE2,'r-d')

title('Consumption')
display('++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
display('chat associated with IRF')
display('++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
%format long;
chat(2) 
% hours plot
subplot(3,3,5)
plot(HRE, 'b-o')
hold on 

plot(zeros(Tend,1))
hold on
plot(Hhat_RE2,'r-d')

title('Hours')

% wage plot
subplot(3,3,6)
plot(wRE, 'b-o')
hold on 
plot(zeros(Tend,1))
hold on
plot(what_RE2,'r-d')
title('Wage')

% rental rate plot
subplot(3,3,7)
plot(RkRE, 'b-o')
hold on 
plot(zeros(Tend,1))
hold on
plot(RKhat_RE2,'r-d')
title('Rental Rate')

% income sum plot
subplot(3,3,8)
plot(SwRE, 'b-o')
hold on 
plot(zeros(Tend,1))

title('Infinite Sum of wage')

% interest rate sum plot
subplot(3,3,9)
plot(SrRE, 'b-o')
hold on 
plot(zeros(Tend,1))

title('Infinite Sum of Interest Rate')

display('xtg')
xtg

figure 

subplot(3,1,1)
set(gca,'fontsize',Font_size)
hp = plot(log(gamma)*ones(1,nsim)*400,'b-o')
set(hp,'linewidth',Line_width)
hold on
plot(y_tg*400, 'r-d')
title('Annualized trend growth beliefs')
legend({'RE','Learning'}, 'Location', 'Southeast','Fontsize',14)

% subplot(2,2,2)
% set(gca,'fontsize',Font_size)
% hpp = plot(yhat*100, 'b-o')  % true output gaps in the learning model
% set(hpp,'linewidth',Line_width)
% hold on
% hpp0 = plot((yhat-xtg)*100, 'r-d')     % subjective estimate of output gaps in the learning model
% set(hpp0,'linewidth',Line_width)
% legend({'True output gaps','Subjective output gaps'}, 'Location', 'Southeast','Fontsize', 13)
% y2hand = get(gca,'ylabel')
% set(y2hand,'string','% dev. from SS','fontsize',17) 
% title('Output gaps: true vs subjective')

subplot(3,1,2)
set(gca,'fontsize',Font_size)
hpp = plot(-xtg*100, 'b-o')  % true output gaps in the learning model
set(hpp,'linewidth',Line_width)
hold on
% hpp0 = plot((yhat-xtg)*100, 'r-d')     % subjective estimate of output gaps in the learning model
% set(hpp0,'linewidth',Line_width)
%legend({'True output gaps','Subjective output gaps'}, 'Location', 'Southeast','Fontsize', 13)
y2hand = get(gca,'ylabel')
set(y2hand,'string','% dev. from SS','fontsize',17) 
title('Bias in output gaps: subjective - true')

subplot(3,1,3)
set(gca,'fontsize',Font_size)
hpp1 = plot(0.53*xtg*100)
set(hpp1,'linewidth',Line_width)
y3hand = get(gca,'ylabel')
set(y3hand,'string','% dev. from SS','fontsize',17) 
% plot()
%  0.53 * 
title('Bias in CAB: subjective - true')

saveas(gcf,'IRF2','epsc')

figure 


subplot(3,1,1)
set(gca,'fontsize',Font_size)
hp2 = plot(-15:nsim, y_trend_o*100, 'b-o')
set(hp2,'linewidth',Line_width)
hold on
plot(-15:nsim, y_trend_s*100, 'r-d')
title('Perceived future and past trend output: two vintages')
legend({'vintage before shock','vintage at the trough'}, 'Location', 'Northwest', 'Fontsize', 10)


subplot(3,1,2)

set(gca,'fontsize',Font_size)

% plot(-15:0,zeros(1,16), 'b-o')  % true output gaps in the learning model
% hold on
plot(-15:0, -xtg_past*100, 'r-d')     % subjective estimate of output gaps in the learning model
%legend({'vintage before the shock', 'vintage at the trough'}, 'Location', 'Southeast','Fontsize', 8)
ythand = get(gca,'ylabel')
set(ythand,'string','% dev. from SS','fontsize',17) 
title('Revision in output gaps before shock: vintage at the trough - vintage before the shock')

subplot(3,1,3)

set(gca,'fontsize',Font_size)

% plot(-15:0,zeros(1,16), 'b-o')  % true output gaps in the learning model
% hold on
plot(-15:0, xtg_past*0.53*100, 'r-d')     % subjective estimate of output gaps in the learning model
%legend({'vintage before the shock', 'vintage at the trough'}, 'Location', 'Southwest','Fontsize', 8)
ythand = get(gca,'ylabel')
set(ythand,'string','% dev. from SS','fontsize',17) 
title('Revision in CAB before shock: vintage at the trough - vintage before the shock')

saveas(gcf,'IRF3','epsc')

figure 

set(gca,'fontsize',Font_size)

plot(-15:0,zeros(1,16), 'b-o')  % true output gaps in the learning model
hold on
plot(-15:0, -xtg_past*100, 'r-d')     % subjective estimate of output gaps in the learning model
legend({'vintage before the shock', 'vintage at the trough'}, 'Location', 'Southeast','Fontsize', 8)
ythand = get(gca,'ylabel')
set(ythand,'string','% dev. from SS','fontsize',17) 
title('Output gaps before shock: two vintages')

figure 

plot(y_tg*400-log(gamma)*ones(1,nsim)*400, 'r-d')
title('Revision in forecast of trend growth beliefs')


%saveas(gcf,'IRF','epsc')
