
%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_excessive(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 ...
coeff_vec_g_RE gy_ratio by_ratio

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

%gainy1 =  0.0045;
% gshock0 = -0.00004;
 %gshock0 = -0.00017;
 
% 2% trend  2009 -- 
% 2008 growth: -1.582215313
cumulative_growth= [-8.103301338
-7.924562056
-8.206075174
-11.0487522
-13.26905276
-13.86740678
-13.80488102
-13.94033763
-13.32663312
];

 % GY ratio 2009 -- 2018
GY_data = 100*(-0.252589845+[0.252589845	0.247929559	0.24089199	0.238873992	0.239054176	0.236577173	0.233198947	0.231290206	0.229903648	0.2299224]);
 % investment data from 2009 to 2018: deviation from trend
inv_data= 100*[-0.015932927
-0.087314077
-0.105038293
-0.110939946
-0.139268419
-0.167703672
-0.178702897
-0.175914113
-0.17636216
-0.171808582
-0.170035855
];

%  euro debt ratio 2009 -- 2018 (relative to steady state 75%)
euro_debtratio = [5.3 10.8	12.7	15.7	17.6	17.8	15.9	15.1	12.7	10.7]; 

% consumption data from 2009 to 2017: deviation from trend: 2007 at 
% steady state
% consumption 2008: 0.009694416
consum_data = 100*[-0.030657383
-0.026281981
-0.022046878
-0.033825518
-0.04851641
-0.055763843
-0.053513553
-0.04976848
-0.037901786
]; 

% working hours : 2009 -- 2014
hours_data = 100*(-99+[99
96
96.4
96
94.1
94.3
])/99;    

% Euro area cab 2009 -- 2015 relative to year 2007; see
% Euro_debtGDPRatio.xlsx
%2017 Spring estimate of 2016; see p. 187 of statistical annex
cab_data = -0.8+1/4*[-3.9	-3.7	-2.6	-1.9	-0.7	-0.2	-0.5 -0.3];

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

shockmat = zeros(1,nsim); 

shockmat(1) =-0.002;  %-0.0024;% std_A;
shockmat(2) = -0.002; %-0.0024;
shockmat(3) = -0.002; %-0.0024;
shockmat(4) = -0.002; %-0.0024;
shockmat(5) = -0.008; %-0.0102;
shockmat(6) = -0.008;%-0.0102;
shockmat(7) = -0.008;%-0.0102;
shockmat(8) = -0.008;%-0.0102;
% shockmat(1) =-0.002;  %-0.0024;% std_A;
% shockmat(2) = -0.002; %-0.0024;
% shockmat(3) = -0.002; %-0.0024;
% shockmat(4) = -0.002; %-0.0024;
% shockmat(5) = -0.008; %-0.0102;
% shockmat(6) = -0.008;%-0.0102;
% shockmat(7) = -0.008;%-0.0102;
% shockmat(8) = -0.008;%-0.0102;


% shockmat(20) =-0.005;  %-0.0024;% std_A;
% shockmat(21) = -0.005; %-0.0024;
% shockmat(22) = -0.005; %-0.0024;
% shockmat(23) = -0.005; %-0.0024;
% shockmat(24) = -0.005; %-0.0102;
% shockmat(25) = -0.005;%-0.0102;
% shockmat(26) = -0.005;%-0.0102;
% shockmat(27) = -0.005;%-0.0102;

 for kkk =1:2
    if kkk==1  %with austerity
        
%         gshock0 = -0.005;%
%         gshock1 = -0.005;
%         gshock2 = -0.005;
%         gshock3 = -0.005;
%         gshock4 = -0.068;
%         gshock5 = -0.068;
%         gshock6 = -0.068;
%         gshock7 = -0.068;
%         gshock8 = -0.08;
%         gshock9 = -0.08;
%         gshock10 = -0.08;
%         gshock11 = -0.08;
        
cum_shock_with = zeros(n1,nsim);

cum_shock_with(1) = shockmat(1);

for ii = 2:nsim
   if ii<= 8 
      cum_shock_with(ii) = cum_shock_with(ii-1) + shockmat(ii);
   else
       cum_shock_with(ii) = cum_shock_with(ii-1);
   end
end


    elseif kkk ==2    %without austerity
        
cum_shock_without = zeros(n1,nsim);

cum_shock_without(1) = shockmat(1);

for ii = 2:nsim
   if ii<= 8 
      cum_shock_without(ii) = cum_shock_without(ii-1) + shockmat(ii);
   else
       cum_shock_without(ii) = cum_shock_without(ii-1);
   end
end

    end
% gshock0 = -0.00;
%  gshock1 = -0.01;
%  gshock2 = -0.01;

% gshock0 = 0;
% gshock1 = 0;
% gshock2 = 0;

Font_size = 11;  
Line_width = 1.5; 



% shockmat(13) = -0.005;
% shockmat(14) = -0.005;
% shockmat(15) = -0.005;

grow_err = zeros(1,nsim); % growth forecast errors

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_ini = eye(3); 
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); 


logG0 = 0;    % initialize government spending
logY0 = log(1/gy_ratio);  %initialize aggregate output


%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
bbhat = zeros(n1,nsim); %budget balance
og_shat = zeros(n1,nsim); %subjective output gaps
cab = zeros(n1,nsim);  % structural balances
og_shat = zeros(n1,nsim);  %

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*inv(mom_mat_ini)*x_vec_m1*(Dlogk0-coeff_mat_k_AL_ini'*x_vec_m1);
           coeff_mat_b_AL(:,t) = coeff_mat_b_AL_ini + gainb*inv(mom_mat_ini)*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*inv(mom_mat_ini)*x_vec_m1*(DlogR0-coeff_mat_R_AL_ini'*x_vec_m1) ;
           
           coeff_mat_w_AL(:,t) = coeff_mat_w_AL_ini + gainw*inv(mom_mat_ini)*x_vec_m1*(Dlogw0-coeff_mat_w_AL_ini'*x_vec_m1) ;
           y_tg(t) = y_tg0; 
           grow_err(t) = 0;  

           
       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*inv(mom_mat_ini)*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*inv(mom_mat_ini)*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*inv(mom_mat_ini)*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*inv(mom_mat_ini)*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)); 
           grow_err(t) = 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* inv(mom_mat_kb(:, 1+3*(t-1):3*t))*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* inv(mom_mat_kb(:, 1+3*(t-1):3*t))*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*inv(mom_mat_kb(:, 1+3*(t-1):3*t))*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*inv(mom_mat_kb(:, 1+3*(t-1):3*t))*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)); 
           grow_err(t) = 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 kkk==1 
      %  if (t==1) | (t==2) | (t==3) | (t==4) | (t==5) | (t==6) | (t==7)  | (t==8) | (t==9)| (t==10)
      if t<=8
              gshock = 0;                     
      else
%          if (cab(t-1)+cab(t-2)+cab(t-3)+cab(t-4))<-0.005 
%              gshock = lamda3*(cab(t-1)+cab(t-2)+cab(t-3)+cab(t-4)+0.005);       

             if (cab(t-5)+cab(t-6)+cab(t-7)+cab(t-8))<-0.005 
          gshock = lamda3*((cab(t-5)+cab(t-6)+cab(t-7)+cab(t-8))+0.005);       
         %     if (cab(t-5)+cab(t-6)+cab(t-7)+cab(t-8))/4<-0.005 
      %   gshock = lamda3*((cab(t-5)+cab(t-6)+cab(t-7)+cab(t-8))/4+0.005); 
           %    if (cab(t-7)+cab(t-8)+cab(t-9)+cab(t-10))/4 < -0.005 m,
          %       gshock = lamda3*((cab(t-7)+cab(t-8)+cab(t-9)+cab(t-10))/4+0.005); 
  
         %gshock = lamda3*(cab(t-5)+0.005); 
        %      gshock = -lamda3*((cab(t-1)+cab(t-2)+cab(t-3)+cab(t-4))/4+0.005)^2*70-0.005; 
%               if t >= 4 && t<= 10 
%               gshock = gshock +0.035;
%               end
%                if t >= 16
%                    gshock = gshock -0.05;
%                end
              else 
                  gshock = 0;
              end
           end
          else 
              gshock = 0; 
          end
           if t==1
              ini_con = [shockmat(t) khat1 bhat1 xtg(t) gshock] ; 
           else
              ini_con = [shockmat(t) khat(t-1) bhat(t-1) xtg(t) gshock]; 
%                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_gshock(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
%             bbhat(t) = -by_ratio * (bhat(t)-yhat(t)-1/gamma*(-shockmat(t)-yhat(t)));
              
              bbhat(t) = -by_ratio *  (1+bhat(t)-1/(gamma*exp(shockmat(t))))/(1+yhat(t));
          else
        %      bbhat(t) = -by_ratio* (bhat(t)-yhat(t)-1/gamma*(bhat(t-1)-shockmat(t)-yhat(t)));
              bbhat(t) = -by_ratio *  (1+bhat(t)-(1+bhat(t-1))/(gamma*exp(shockmat(t))))/(1+yhat(t));
          end
          
           og_shat(t) = yhat(t) - xtg(t);
           cab(t) = bbhat(t) - 0.53* og_shat(t); % structural balance ratio
           
          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); 
           DlogG(t) = ghat(t)-0+log(gamma)+shockmat(t);  % added government expenditures
           logG(t) = DlogG(t) + logG0; 
           logY(t) = DlogY(t) + logY0; 
          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);
           DlogG(t) = ghat(t)-ghat(t-1)+log(gamma)+shockmat(t); % added government expenditures
           logG(t) = DlogG(t) + logG(t-1);  
           logY(t) = DlogY(t) + logY(t-1); 
           
          end

    

    end % end Tend loop over simul periods
   
 %    [y_tg_min , index] = min(y_tg);
%      display('index')
%      index

    Yg2008 = sum(DlogY(1:4));
    Yg2009 = sum(DlogY(5:8));
    display('GDP growth 2008')
    Yg2008
    display('GDP growth 2009')
    Yg2009
    

% display('output growth in 2008')
% (exp(logY(4))+exp(logY(3))+exp(logY(2))+exp(logY(1)))/(exp(logY0)+exp(logY0-log(gamma))+exp(logY0-2*log(gamma))+exp(logY0-3*log(gamma))) - 1
% 
% display('output growth in 2009')
% (exp(logY(8))+exp(logY(7))+exp(logY(6))+exp(logY(5)))/(exp(logY(4))+exp(logY(3))+exp(logY(2))+exp(logY(1))) - 1
  index = 13; 
    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

    
       %%% change in government spending GDP ratios added 2 August 2019
       G2009 = exp(logG(5))+exp(logG(6))+exp(logG(7))+exp(logG(8)); 
       Y2009 = exp(logY(5))+exp(logY(6))+exp(logY(7))+exp(logY(8)); 
       GY2009 = G2009/Y2009;      
       
       G2010 = exp(logG(9))+exp(logG(10))+exp(logG(11))+exp(logG(12)); 
       Y2010 = exp(logY(9))+exp(logY(10))+exp(logY(11))+exp(logY(12)); 
       GY2010 = G2010/Y2010; 
       
       G2011 = exp(logG(13))+exp(logG(14))+exp(logG(15))+exp(logG(16)); 
       Y2011 = exp(logY(13))+exp(logY(14))+exp(logY(15))+exp(logY(16)); 
       GY2011 = G2011/Y2011; 
       
       G2012 = exp(logG(17))+exp(logG(18))+exp(logG(19))+exp(logG(20)); 
       Y2012 = exp(logY(17))+exp(logY(18))+exp(logY(19))+exp(logY(20)); 
       GY2012 = G2012/Y2012; 

       display('GY2009')
       GY2009
       
       display('GY2010')
       GY2010
       
       display('GY2011')
       GY2011
       
       display('GY2012')
       GY2012

       
       
       display('GY2010-GY2009')
       GY2010-GY2009
       
       display('GY2011-GY2009')
       GY2011-GY2009
       
       display('GY2012-GY2009')
       GY2012-GY2009

       
       
% 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);
bbRE = 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
%     display('time period')
%     j
                         
%            elseif t==4
%               gshock = gshock0;                
%            elseif t==6
%               lam2shock = lam2shock0; 
%            elseif t==8
%               gshock = gshock1;  
%            elseif t==9
%               gshock = gshock2;  
%            elseif t==10
%               gshock = gshock2;  
%            elseif t==11
%               gshock = gshock2;                
        
              gshock_re = 0; 
        
       if j==1    
            HRE(j) = coeff_vec_s_RE(5)*shockmat(j) + coeff_vec_g_RE(5)*gy_ratio*gshock_re;
            yyRE(j) = coeff_vec_s_RE(3)*shockmat(j)  + coeff_vec_g_RE(3)*gy_ratio*gshock_re; 
            cRE(j) =coeff_vec_s_RE(6)*shockmat(j) + coeff_vec_g_RE(6)*gy_ratio*gshock_re;
            wRE(j) =coeff_vec_s_RE(2)*shockmat(j) + coeff_vec_g_RE(2)*gy_ratio*gshock_re;
            caputRE(j) = coeff_vec_s_RE(8)*shockmat(j) + coeff_vec_g_RE(8)*gy_ratio*gshock_re; 
      %      RtRE(j,i) = coeff_vec_s_RE(4)*shockmat(j,i);
            REinv(j) = coeff_vec_s_RE(7)*shockmat(j) + coeff_vec_g_RE(7)*gy_ratio*gshock_re;
            kRE(j) = coeff_vec_s_RE(9)*shockmat(j) + coeff_vec_g_RE(9)*gy_ratio*gshock_re;
            gtRE(j) = coeff_vec_s_RE(4)*shockmat(j) + coeff_vec_g_RE(4)*gy_ratio*gshock_re;
      %      atRE(j,i) = coeff_vec_s_RE(12)*shockmat(j,i);
            btRE(j) = coeff_vec_s_RE(11)*shockmat(j) + coeff_vec_g_RE(11)*gy_ratio*gshock_re;
            RkRE(j) = coeff_vec_s_RE(1)*shockmat(j) + coeff_vec_g_RE(1)*gy_ratio*gshock_re;

                ini_con_RE = [shockmat(j) 0  0 gshock_re]; 

               [y_vec1_RE, fore_dis_RE]  = Learningsol_RE_gshock(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) + coeff_vec_g_RE(5)*gy_ratio*gshock_re;
            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)+ coeff_vec_g_RE(3)*gy_ratio*gshock_re;          
            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)+ coeff_vec_g_RE(6)*gy_ratio*gshock_re;         
            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)+ coeff_vec_g_RE(2)*gy_ratio*gshock_re;                      
            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)+ coeff_vec_g_RE(7)*gy_ratio*gshock_re;                      
            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)+ coeff_vec_g_RE(9)*gy_ratio*gshock_re; 
            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)+ coeff_vec_g_RE(4)*gy_ratio*gshock_re;
         %   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)+ coeff_vec_g_RE(11)*gy_ratio*gshock_re;
            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)+ coeff_vec_g_RE(1)*gy_ratio*gshock_re;
   
            
               ini_con_RE = [shockmat(j) kRE(j-1) btRE(j-1) gshock_re]; 
%                display('ini_con_RE')
%                ini_con_RE
               [y_vec1_RE, fore_dis_RE]  = Learningsol_RE_gshock(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);
%             if t==1
%               bbhat(t) = -by_ratio* (bhat(t)-yhat(t)-1/gamma*(-shockmat(t)-yhat(t)));
%           else
%               bbhat(t) = -by_ratio* (bhat(t)-yhat(t)-1/gamma*(bhat(t-1)-shockmat(t)-yhat(t)));
%           end        
%            

%           if t==1
% %             bbhat(t) = -by_ratio * (bhat(t)-yhat(t)-1/gamma*(-shockmat(t)-yhat(t)));
%               
%               bbhat(t) = -by_ratio *  (1+bhat(t)-1/(gamma*exp(shockmat(t))))/(1+yhat(t));
%           else
%         %      bbhat(t) = -by_ratio* (bhat(t)-yhat(t)-1/gamma*(bhat(t-1)-shockmat(t)-yhat(t)));
%               bbhat(t) = -by_ratio *  (1+bhat(t)-(1+bhat(t-1))/(gamma*exp(shockmat(t))))/(1+yhat(t));
%           end
%           
          
          if j==1
              bbRE(j) =   -by_ratio* (1+btRE(j)-1/(gamma*exp(shockmat(j))))/(1+yyRE(j));
              
       %       bbRE(j) =   -by_ratio* (1+btRE(j)-yyRE(j)-1/gamma*(-shockmat(j)-yyRE(j)))  ;
          else
     %         bbRE(j) =  -by_ratio* (btRE(j)-yyRE(j)-1/gamma*(btRE(j-1)-shockmat(t)-yyRE(j))) ;
                     bbRE(j) =   -by_ratio* (1+btRE(j)-(1+btRE(j-1))/(gamma*exp(shockmat(j))))/(1+yyRE(j));
          end
   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')

if kkk ==1 
   y_RE_with = yyRE*100+cum_shock_with*100;
   y_with = yhat*100+cum_shock_with*100;
   inv_RE_with = REinv*100+cum_shock_with*100;
   inv_with = ihat*100+cum_shock_with*100;
   c_RE_with = cRE*100+cum_shock_with*100;
   c_with = chat*100+cum_shock_with*100;
   H_RE_with =HRE*100;
   H_with = Hhat*100;
   g_RE_with = gtRE*100+cum_shock_with*100;
   g_RE_with1 =  gy_ratio * (gtRE - yyRE)*100; % government expenditure GDP ratio under RE with austerity
   g_with = ghat*100+cum_shock_with*100;    
   g_with1 = gy_ratio * (ghat - yhat)*100; % government expenditures GDP ratio under learning with austerity
   by_RE_with = (btRE-yyRE)*100;
   by_with = (bhat-yhat)*100;
   cab_with = cab;
   
 % used to calculate spending multiplier
   yhat_A = yhat;  % output under learning and with austerity 
   ghat_A = ghat;  % government spending under learning and with austerity
   yhat_A_RE = yyRE;  % output under RE and with austerity
   ghat_A_RE = gtRE;       % government spending under RE and with austerity
   
elseif kkk==2
   y_RE_without = yyRE*100+cum_shock_without*100;
   y_without = yhat*100+cum_shock_without*100;
   inv_RE_without = REinv*100+cum_shock_without*100;
   inv_without = ihat*100+cum_shock_without*100;
   c_RE_without = cRE*100+cum_shock_without*100;
   c_without = chat*100+cum_shock_without*100;
   H_RE_without =HRE*100;
   H_without = Hhat*100;
   g_RE_without = gtRE*100+cum_shock_without*100;
   g_RE_without1 = gy_ratio * (gtRE-yyRE)*100;    % government expenditure GDP ratio under RE without austerity
   g_without = ghat*100+cum_shock_without*100;  % 
   g_without1 = gy_ratio * (ghat - yhat) * 100; % government expenditures GDP ratio under learning without austerity
   by_RE_without = (btRE-yyRE)*100;
   by_without = (bhat-yhat)*100;
   cab_without =cab;
   
% used to calculate spending multiplier
   
  yhat_noA = yhat;
  ghat_noA = ghat; 
  yhat_noA_RE = yyRE;
  ghat_noA_RE = gtRE; 
  
% display('yhat_A(9:NNN) - yhat_noA(9:NNN)')
% yhat_A(9:40) - yhat_noA(9:40)
% 
% display('ghat_A(9:NNN) - ghat_noA(9:NNN)')
% ghat_A(9:40) - ghat_noA(9:40)
% 
% 
% display('yhat_A_RE(9:NNN) - yhat_noA_RE(9:NNN)')
% yhat_A_RE(9:40) - yhat_noA_RE(9:40)
% 
% display('ghat_A_RE(9:NNN) - ghat_noA_RE(9:NNN)')
% ghat_A_RE(9:40) - ghat_noA_RE(9:40)


  for NNN = 9: nsim
  
  multiplier(NNN) = sum(yhat_A(9:NNN) - yhat_noA(9:NNN))/sum(ghat_A(9:NNN) - ghat_noA(9:NNN))/gy_ratio; % multiplier under leanring 
  multiplier_RE(NNN) = sum(yhat_A_RE(9:NNN) - yhat_noA_RE(9:NNN))/sum(ghat_A_RE(9:NNN) - ghat_noA_RE(9:NNN))/gy_ratio;  % multiplier under RE 
  
  multiplier1(NNN) = sum(yhat_A(9:NNN))/sum(ghat_A(9:NNN))/gy_ratio;  
  multiplier_RE1(NNN) = sum(yhat_A_RE(9:NNN))/sum(ghat_A_RE(9:NNN))/gy_ratio;
  end
  display('+++++++++++++++++++multiplier++++++++++++++++')
  
  multiplier
  multiplier1
  multiplier_RE
  multiplier_RE1
  
%    
% set(gca,'fontsize',Font_size)
% hi1 = plot(g_RE_with1, 'b-o')
% set(hi1,'linewidth',Line_width)  
% 
% hold on 
% plot(g_with1,'r-d')
% hold on 
% hi1 = plot(g_RE_without1, 'b-*')
% set(hi1,'linewidth',Line_width)
% 
% hold on 
% plot(g_without1,'r-x')
% hold on 
% plot(zeros(Tend,1))
%legend({'RE', 'Learning'}, 'Location', 'Southeast', 'FontSize', 15);   
end


if kkk==2
%     figure 
% subplot(3,1,1)
% plot(bbRE_with*100, 'b-o')
% hold on
% plot(bbhat_with*100, 'r-d')
% title('Budget balance GDP ratio')
% plot(bbRE_without*100, 'b-o')
% hold on
% plot(bbhat_without*100, 'r-d')
% title('Budget balance GDP ratio')
% 
% hold on 
% plot(zeros(Tend,1))
% legend({'RE -austerity', 'Learning - austerity','RE', 'Learning'}, 'Location', 'Southeast', 'FontSize', 9);

% subplot(3,1,2)
% plot(yyRE_with*100, 'b-o')
% hold on 
% plot(og_shat_with*100, 'r-d')
% hold on 
% plot(yyRE_without*100, 'b-o')
% hold on 
% plot(og_shat_without*100, 'r-d')
% hold on 
% 
% plot(zeros(Tend,1))
% legend({'RE -austerity', 'Learning - austerity','RE', 'Learning'}, 'Location', 'Southeast', 'FontSize', 9);
% title('Output gaps')


% cumulative_growth= [ -0.582215313
% -6.103301338
% -4.924562056
% -4.206075174
% -6.048752201
% -7.269052757
% -6.86740678
% -5.804881015
% -4.940337627
% -3.326633118
% ];
% cumulative_growth= [-1.082215313
% -7.103301338
% -6.424562056
% -6.206075174
% -8.548752201
% -10.26905276
% -10.36740678
% -9.804881015
% -9.440337627
% -8.326633118
% -7.994006527
% 
% 
% ];

years = 2009:2013
numyears = length(years);

jan1s = zeros(numyears,1);
quarters = zeros(4,numyears);
for Yearno = 1:numyears
 thisyear = years(Yearno);
 jan1s(Yearno) = datenum([thisyear 1 1]);
 quarters(1,Yearno) = jan1s(Yearno);
 quarters(2,Yearno) = datenum([thisyear 4 1]);
 quarters(3,Yearno) = datenum([thisyear 7 1]);
 quarters(4,Yearno) = datenum([thisyear 10 1]);
end

figure
subplot(2,1,1)
set(gca,'fontsize',Font_size)
h1 = plot(quarters(:),RKhat(5:Tend),'r-')
set(h1,'linewidth',Line_width)
set(gca,'XTick', jan1s);
datetick('x', 'yyyy QQ');
grid on
grid minor
hold on 
title('interest rate')
subplot(2,1,2)
set(gca,'fontsize',Font_size)
h1 = plot(quarters(:),what(5:Tend),'r-')
set(h1,'linewidth',Line_width)
set(gca,'XTick', jan1s);
datetick('x', 'yyyy');
grid on
grid minor
hold on 
title('wage rate')


figure 
subplot(3,2,1)
set(gca,'fontsize',Font_size)
h1 = plot(quarters(linspace(4,Tend-4,Tend/4-1)),cumulative_growth(1:Tend/4-1),'g-')
set(h1,'linewidth',Line_width)

set(gca,'XTick', jan1s);
datetick('x', 'yyyy');
hold on 
h1 = plot(quarters(:),y_RE_with(5:Tend), 'b-')
set(h1,'linewidth',Line_width)

grid on
grid minor
hold on 
plot(quarters(:),y_with(5:Tend),'r-')
y_with

% hold on 
% plot(quarters(:),y_without,'r--')

hold on 
plot(quarters(:),zeros(Tend-4,1))

ylhand = get(gca,'ylabel')
set(ylhand,'string','% dev. from unshocked BGP','fontsize',10) 
hold on 
legend({'data','RE', 'Learning'}, 'Location', 'Southwest', 'FontSize', 4.5);
legend boxoff  

title('Output')
% display('average output loss over three years')
% mean(y_with(9:20))- mean(y_without(9:20))

% display('cumulative output loss over six years')
% (mean(y_with(9:nsim-8))- mean(y_without(9:nsim-8)))*6



% capital stock plot

%plot(zeros(Tend,1))

%output plot

% investment plot
subplot(3,2,2)
set(gca,'fontsize',Font_size)
h1 = plot(quarters(linspace(4,Tend-4,Tend/4-1)),inv_data(1:Tend/4-1)-inv_data(1)+inv_with(8),'g-')
set(h1,'linewidth',Line_width)

set(gca,'XTick', jan1s);
datetick('x', 'yyyy');
hold on 

h2 = plot(quarters(:),inv_RE_with(5:Tend), 'b-')
set(h2,'linewidth',Line_width)
set(gca,'XTick', jan1s);
datetick('x', 'yyyy');
grid on
grid minor
hold on 

plot(quarters(:),inv_with(5:Tend),'r-')
hold on
% h3 = plot(quarters(:),inv_RE_without(5:Tend), 'b--')
% set(h3,'linewidth',Line_width)
% % plot(quarters(:),inv_without,'r--')
% hold on
plot(zeros(Tend,1))

title('Investment')

% consumption plot
subplot(3,2,3)
set(gca,'fontsize',Font_size)
h1 = plot(quarters(linspace(4,Tend-4,Tend/4-1)),consum_data(1:Tend/4-1)-consum_data(1)+c_with(8),'g-')
set(h1,'linewidth',Line_width)

set(gca,'XTick', jan1s);
datetick('x', 'yyyy');
hold on 
h3 = plot(quarters(:),c_RE_with(5:Tend), 'b-')
set(h3,'linewidth',Line_width)
set(gca,'XTick', jan1s);
datetick('x', 'yyyy');
grid on
grid minor
hold on 
plot(quarters(:),c_with(5:Tend),'r-')
hold on 
% h4 = plot(quarters(:),c_RE_without, 'b--')
%set(h4,'linewidth',Line_width)
hold on 
% plot(quarters(:),c_without,'r--')
% hold on 
plot(quarters(:),zeros(Tend-4,1))

title('Consumption')
% display('++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
% display('chat associated with IRF')
% display('++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
% %format long;
% chat(2) 
% hours plot
subplot(3,2,4)

set(gca,'fontsize',Font_size)
h1 = plot(quarters(linspace(4,Tend-4,Tend/4-1)),hours_data(1:Tend/4-1)-hours_data(1)+H_with(8),'g-')
set(h1,'linewidth',Line_width)

set(gca,'XTick', jan1s);
datetick('x', 'yyyy');
hold on 

h4 = plot(quarters(:),H_RE_with(5:Tend), 'b-')
set(gca,'XTick', jan1s);
datetick('x', 'yyyy');
grid on
grid minor
set(h4,'linewidth',Line_width)
hold on 
plot(quarters(:),H_with(5:Tend),'r-')
% hold on 
% h5 = plot(quarters(:),H_RE_without, 'b--')
% set(h5,'linewidth',Line_width)
hold on 
% plot(quarters(:),H_without,'r--')
% hold on 

plot(quarters(:),zeros(Tend-4,1))

title('Hours')


%GY_data= 0.5+100*[-0.004660286	-0.011697856	-0.013715853	-0.013535669	-0.016012672	-0.019390898];




% government expenditures
subplot(3,2,5)
set(gca,'fontsize',Font_size)
h1 = plot(quarters(linspace(4,Tend-4,Tend/4-1)),GY_data(1:Tend/4-1)-GY_data(1)+g_with1(8),'g-')
set(h1,'linewidth',Line_width)

set(gca,'XTick', jan1s);
datetick('x', 'yyyy');
hold on 

hi1 = plot(quarters(:),g_RE_with1(5:Tend), 'b-')
set(gca,'XTick', jan1s);
datetick('x', 'yyyy');
grid on
grid minor
set(hi1,'linewidth',Line_width)  

hold on 
plot(quarters(:),g_with1(5:Tend),'r-')
hold on 
% hi1 = plot(quarters(:),g_RE_without1(5:Tend), 'b--')
% set(hi1,'linewidth',Line_width)
% h1=plot(quarters(linspace(12,Tend, Tend/4-2)), GY_data(1:end-2),'g-')
% set(h1,'linewidth',Line_width)
% hold on 
% plot(quarters(:),g_without1,'r--')
% hold on 
plot(quarters(:),zeros(Tend-4,1))
%legend({'RE', 'Learning'}, 'Location', 'Southeast', 'FontSize', 15);

title('Gov. exp. ratio')
display('g_with1')

g_with1


% subplot(3,2,5)
% set(gca,'fontsize',Font_size)
% hi1 = plot(gy_ratio*(gtRE - yyRE) * 100, 'b-o')
% set(hi1,'linewidth',Line_width)
% 
% hold on 
% plot(gy_ratio*(ghat- yhat)*100,'r-d')
% hold on 
% plot(zeros(Tend,1))
% %legend({'RE', 'Learning'}, 'Location', 'Southeast', 'FontSize', 15);
% 
% title('Gov. exp. ratio')

%plot(zeros(Tend,1))

%output plot

%  euro debt ratio 2010 -- 2018 (relative to 2009)
% euro_debtratio = [5.5	7.4	10.4	12.3	12.5	10.6	9.8	7.4	5.4];




% euro debt ratio 2008 -- (relative to 2007)
%euro_debtratio1 = [3.7	14.3	19.8	21.7	24.7	26.6	26.8	24.9	24.1	21.7	19.7];

subplot(3,2,6)
set(gca,'fontsize',Font_size)

h1 = plot(quarters(linspace(4,Tend-4,Tend/4-1)),euro_debtratio(1:Tend/4-1)-euro_debtratio(1)+by_with(8),'g-')
set(h1,'linewidth',Line_width)

set(gca,'XTick', jan1s);
datetick('x', 'yyyy');
hold on 

h2 = plot(quarters(:),by_RE_with(5:Tend), 'b-')
set(gca,'XTick', jan1s);
datetick('x', 'yyyy');
grid on
grid minor
set(h2,'linewidth',Line_width)
hold on 
plot(quarters(:),by_with(5:Tend),'r-')
hold on
% h2 = plot(quarters(:),by_RE_without, 'b--')
% set(h2,'linewidth',Line_width)
% hold on 
plot(quarters(:),zeros(Tend-4,1))
hold on 
% plot(quarters(:),by_without,'r--')
% hold on 
% % plot(quarters(:),zeros(Tend-4,1))
% % hold on 
% h1 = plot(quarters(linspace(8,(Tend-4),(Tend/4-2))),euro_debtratio(1:(Tend/4-2)),'g-')
% set(h1,'linewidth',Line_width)
% h1 = plot(quarters(linspace(4,(Tend-4),(Tend/4-1))),euro_debtratio1(1:(Tend/4-1)),'g-')
% set(h1,'linewidth',Line_width)
title('Debt/GDP ratio')
saveas(gcf,'IRF_excessive','epsc')
end

if kkk ==1 
    trg_RE = log(gamma)*ones(1,nsim)*400;
    trg_with = y_tg*400;  %adaptive learning with austerity
    og_bias_with =-xtg*100; 
    cab_bias_with =0.53*xtg*100; 
elseif kkk == 2
    trg_without = y_tg*400; % without austerity
    og_bias_without =-xtg*100; 
    cab_bias_without =0.53*xtg*100; 
end

if kkk==2
    
figure 

subplot(2,2,1)
set(gca,'fontsize',Font_size)
hp = plot(quarters(:),trg_RE(5:Tend),'b-o')
set(gca,'XTick', jan1s);
datetick('x', 'yyyy');
set(hp,'linewidth',Line_width)
grid on
grid minor
hold on
plot(quarters(:),trg_with(5:Tend), 'r-d')
% hold on
% plot(quarters(:),trg_without, 'r-x')
legend({'RE','Learning'}, 'Location', 'Southeast','Fontsize',4)
legend boxoff  
% h = findobj('type', 'axes');  % Find all sets of axes
% set(h(1), 'visible', 'off')    % Hides the legend's axes (legend border and background)
title('Annualized trend growth beliefs')


% 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(2,2,2)
set(gca,'fontsize',Font_size)
hpp = plot(quarters(:),og_bias_with(5:Tend), 'r-d')  % true output gaps in the learning model
set(hpp,'linewidth',Line_width)
set(gca,'XTick', jan1s);
datetick('x', 'yyyy');
% hold on
% plot(quarters(:),og_bias_without, 'r-*')  % 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','percentage','fontsize',12) 
grid on
grid minor
title('Estimation error in OG')
%: subjective - true
subplot(2,2,3)
set(gca,'fontsize',Font_size)
hpp1 = plot(quarters(:),cab_bias_with(5:Tend), 'r-d')
set(hpp1,'linewidth',Line_width)
set(gca,'XTick', jan1s);
datetick('x', 'yyyy');
% hold on
% plot(quarters(:),cab_bias_without, 'r-*')
y3hand = get(gca,'ylabel')
set(y3hand,'string','percentage','fontsize',12) 
grid on
grid minor
% plot()
%  0.53 * 
title('Estimation error in CAB')
%: subjective - true


subplot(2,2,4)
% hpp0 = plot(zeros(n1,nsim), 'b-o')
% set(hpp0,'linewidth',Line_width)
% hold on


%h1 = plot(quarters(linspace(4,Tend-4,Tend/4-1)),cab_data(1:Tend/4-1)-cab_data(1)+100*cab_with(8),'g-')
h1 = plot(quarters(linspace(4,Tend-4,Tend/4-1)),cab_data(1:Tend/4-1),'g-')
set(h1,'linewidth',Line_width)

set(gca,'XTick', jan1s);
datetick('x', 'yyyy');
hold on 

% h1 = plot(quarters(linspace(4,Tend,Tend/4)),cab_data(1:(Tend/4)),'g-')
% set(h1,'linewidth',Line_width)
% hold on
% set(gca,'XTick', jan1s);
% datetick('x', 'yyyy');


plot(quarters(:),cab_with(5:Tend)*100, 'r-d')
% hold on 
% plot(quarters(:),cab_without*100, 'r-x')
grid on
grid minor
legend({'data','Learning'}, 'Location', 'North','Fontsize',4)
legend boxoff  
plot(quarters(:),zeros(Tend-4,1))
hold on
%hold on
%legend({'Learning - austerity', 'Learning','RE'}, 'Location', 'Southeast', 'FontSize', 8);

title('CAB')


saveas(gcf,'IRF2_excessive','epsc')
end
figure 


years1 = 2004:2007
numyears1 = length(years1);

jan1s1 = zeros(numyears1,1);
quarters1 = zeros(4,numyears1);
for Yearno1 = 1:numyears1
 thisyear1 = years1(Yearno1);
 jan1s1(Yearno1) = datenum([thisyear1 1 1]);
 quarters1(1,Yearno1) = jan1s1(Yearno1);
 quarters1(2,Yearno1) = datenum([thisyear1 4 1]);
 quarters1(3,Yearno1) = datenum([thisyear1 7 1]);
 quarters1(4,Yearno1) = datenum([thisyear1 10 1]);
end
quarters1

years2 = 2004:2012
numyears2 = length(years2);

jan1s2 = zeros(numyears2,1);
quarters2 = zeros(4,numyears2);
for Yearno2 = 1:numyears2
 thisyear2 = years2(Yearno2);
 jan1s2(Yearno2) = datenum([thisyear2 1 1]);
 quarters2(1,Yearno2) = jan1s2(Yearno2);
 quarters2(2,Yearno2) = datenum([thisyear2 4 1]);
 quarters2(3,Yearno2) = datenum([thisyear2 7 1]);
 quarters2(4,Yearno2) = datenum([thisyear2 10 1]);
end
quarters2


subplot(3,1,1)
set(gca,'fontsize',Font_size)
hp2 = plot(quarters2(:), y_trend_o(1:36)*100, 'r-o')
set(gca,'XTick', jan1s2);
datetick('x', 'yyyy');
set(hp2,'linewidth',Line_width)
hold on
grid on
grid minor
plot(quarters2(:), y_trend_s(1:36)*100, 'r-d')
title('Perceived future and past trend output: two vintages')
legend({'pre-crisis vintage (08 Spring)','vintage 10 Spring'}, 'Location', 'Northwest', 'Fontsize', 8)


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(quarters1(:), -xtg_past*100, 'r-d')     % subjective estimate of output gaps in the learning model
set(gca,'XTick', jan1s1);
datetick('x', 'yyyy');
%legend({'vintage before the shock', 'vintage at the trough'}, 'Location', 'Southeast','Fontsize', 8)
ythand = get(gca,'ylabel')
set(ythand,'string','%','fontsize',15) 
grid on
grid minor
title('Revision to output gaps before shock: vintage 2010 Spring - vintage 2008 Spring')

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(quarters1(:), xtg_past*0.53*100, 'r-d')     % subjective estimate of output gaps in the learning model
set(gca,'XTick', jan1s1);
datetick('x', 'yyyy');
%legend({'vintage before the shock', 'vintage at the trough'}, 'Location', 'Southwest','Fontsize', 8)
ythand = get(gca,'ylabel')
set(ythand,'string','%','fontsize',15) 
grid on
grid minor
title('Revision to CAB before shock: vintage 2010 Spring - vintage 2008 Spring')

saveas(gcf,'IRF3_excessive','epsc')

figure 

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

display('change in output gaps')
yhat(6)-yhat(5)
display('change in debt GDP ratio')
(bhat(6)-yhat(6))- (bhat(5)-yhat(5)) 

display('grow_err(6)')
grow_err(6)

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