
%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;
gainy = 0.003; 
semi_ela = 0.55; % semi-elasticity
lamda3 = 1.1;
cumulative_growth= [-2.250410105
-8.932462475
-8.445939465
-9.069311766
-13.08833353
-16.03652159
-17.57968887
-18.02015283
-18.12970937
];
GY_data = 0*[-0.250720991	-0.634253978	-0.61275509	-0.535323485	-0.671539461	-1.003472257];
euro_debtratio = 0* [10.8	12.7	15.7	17.6	17.8	15.9	15.1	12.7	10.7]; 
cab_data = 0*[-2.2	-3.9	-3.7	-2.6	-1.9	-0.7	-0.2	-0.5 -0.3];
 

 
%% GENERATES SHOCKS
n1 =1;
 nsim =40;
  Tend = nsim;

shockmat = zeros(1,nsim); 


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

 n1 = 1;
 growth_q1_f = zeros(n1, nsim); % 1Q-ahead growth forecast
growth_q2_f = zeros(n1, nsim); % 2Q-ahead growth forecast
growth_q3_f = zeros(n1, nsim); % 3Q-ahead growth forecast
growth_q4_f = zeros(n1, nsim); % 4Q-ahead growth forecast
growth_q5_f = zeros(n1, nsim); % 5Q-ahead growth forecast
growth_q6_f = zeros(n1, nsim); % 6Q-ahead growth forecast
growth_q7_f = zeros(n1, nsim); % 7Q-ahead growth forecast
growth_q8_f = zeros(n1, nsim); % 8Q-ahead growth forecast
growth_q9_f = zeros(n1, nsim); % 9Q-ahead growth forecast
growth_q10_f = zeros(n1, nsim); % 10Q-ahead growth forecast
growth_q11_f = zeros(n1, nsim); % 11Q-ahead growth forecast
growth_q12_f = zeros(n1, nsim); % 12Q-ahead growth forecast


growth_q2_f_k = zeros(2,1);
growth_q3_f_k = zeros(2,1);
growth_q4_f_k = zeros(2,1);

growth_q5_f_k = zeros(2,1);
growth_q6_f_k = zeros(2,1);
growth_q7_f_k = zeros(2,1);
growth_q8_f_k = zeros(2,1);
growth_q9_f_k = zeros(2,1);
growth_q10_f_k = zeros(2,1);
growth_q11_f_k = zeros(2,1);
growth_q12_f_k = zeros(2,1);


 for kkk =1:2
    if kkk==1
        
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
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

Font_size = 11;  
Line_width = 1.5; 



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);
coeff_mat_y_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

Dlogy0 = log(gamma);
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)];
%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 = zeros(n1,nsim);  % belief vector about trend growth

RKhat = zeros(n1,nsim);  %  rental rate
what = zeros(n1,nsim);   % wages 
yhat = zeros(n1,nsim);   %  output
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
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 
           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) ;
           coeff_mat_y_AL(:,t) = coeff_mat_y_AL_ini + gainy*inv(mom_mat_ini)*x_vec_m1*(Dlogy0-coeff_mat_y_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;
           
           
           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);
           coeff_mat_y_AL(:,t) = coeff_mat_y_AL(:,t-1) + gainy*inv(mom_mat_ini)*x_vec_ini*(DlogY(t-1)-coeff_mat_y_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)));
           
           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);
           coeff_mat_y_AL(:,t) = coeff_mat_y_AL(:,t-1) + gainy*inv(mom_mat_kb(:, 1+3*(t-1):3*t))*x_vec*(DlogY(t-1)-coeff_mat_y_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<=5
              gshock = 0;                     
      elseif t == 6
         if cab(t-5)<-0.005 
            gshock = lamda3*(cab(t-5)+0.005);       
         else 
            gshock = 0;
         end

      elseif t==7
             if (cab(t-5)+cab(t-6))<-0.005 
                gshock = lamda3*((cab(t-5)+cab(t-6))+0.005);       
             else 
                  gshock = 0;
             end
         
      elseif t==8 
             if (cab(t-5)+cab(t-6)+cab(t-7))<-0.005 
                gshock = lamda3*((cab(t-5)+cab(t-6)+cab(t-7))+0.005);       
             else 
                gshock = 0;
             end

      else
          
             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);       
             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]; 
           end
           %      xw : trend belief in w equation
           %      xR : trend belief in R equation
           % ini_con: ordering gam lagk
           [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
           
           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 *  (1+bhat(t)-1/(gamma*exp(shockmat(t))))/(1+yhat(t));
          else
              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) - semi_ela* 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);
           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);
           logG(t) = DlogG(t) + logG(t-1);  
           logY(t) = DlogY(t) + logY(t-1);            
          end
       if kkk ==1
                   growth_q2_f_k = [coeff_mat_k_AL(1,t) ;coeff_mat_b_AL(1,t)]+...
                   [coeff_mat_k_AL(2,t) coeff_mat_k_AL(3,t);coeff_mat_b_AL(2,t) coeff_mat_b_AL(3,t)]*[Dlogk(1,t); Dlogb(1,t)];
               
               growth_q3_f_k = [coeff_mat_k_AL(1,t) ; coeff_mat_b_AL(1,t)]+...
                   [coeff_mat_k_AL(2,t) coeff_mat_k_AL(3,t);coeff_mat_b_AL(2,t) coeff_mat_b_AL(3,t)]* growth_q2_f_k;
               growth_q4_f_k = [coeff_mat_k_AL(1,t) ; coeff_mat_b_AL(1,t)]+...
                   [coeff_mat_k_AL(2,t) coeff_mat_k_AL(3,t);coeff_mat_b_AL(2,t) coeff_mat_b_AL(3,t)]* growth_q3_f_k;
               growth_q5_f_k = [coeff_mat_k_AL(1,t) ; coeff_mat_b_AL(1,t)]+...
                   [coeff_mat_k_AL(2,t) coeff_mat_k_AL(3,t);coeff_mat_b_AL(2,t) coeff_mat_b_AL(3,t)]* growth_q4_f_k;

                              growth_q6_f_k = [coeff_mat_k_AL(1,t) ; coeff_mat_b_AL(1,t)]+...
                   [coeff_mat_k_AL(2,t) coeff_mat_k_AL(3,t);coeff_mat_b_AL(2,t) coeff_mat_b_AL(3,t)]* growth_q5_f_k;

                              growth_q7_f_k = [coeff_mat_k_AL(1,t) ; coeff_mat_b_AL(1,t)]+...
                   [coeff_mat_k_AL(2,t) coeff_mat_k_AL(3,t);coeff_mat_b_AL(2,t) coeff_mat_b_AL(3,t)]* growth_q6_f_k;

                              growth_q8_f_k = [coeff_mat_k_AL(1,t) ; coeff_mat_b_AL(1,t)]+...
                   [coeff_mat_k_AL(2,t) coeff_mat_k_AL(3,t);coeff_mat_b_AL(2,t) coeff_mat_b_AL(3,t)]* growth_q7_f_k;

                              growth_q9_f_k = [coeff_mat_k_AL(1,t) ; coeff_mat_b_AL(1,t)]+...
                   [coeff_mat_k_AL(2,t) coeff_mat_k_AL(3,t);coeff_mat_b_AL(2,t) coeff_mat_b_AL(3,t)]* growth_q8_f_k;

                              growth_q10_f_k = [coeff_mat_k_AL(1,t) ; coeff_mat_b_AL(1,t)]+...
                   [coeff_mat_k_AL(2,t) coeff_mat_k_AL(3,t);coeff_mat_b_AL(2,t) coeff_mat_b_AL(3,t)]* growth_q9_f_k;

                              growth_q11_f_k = [coeff_mat_k_AL(1,t) ; coeff_mat_b_AL(1,t)]+...
                   [coeff_mat_k_AL(2,t) coeff_mat_k_AL(3,t);coeff_mat_b_AL(2,t) coeff_mat_b_AL(3,t)]* growth_q10_f_k;

                              growth_q12_f_k = [coeff_mat_k_AL(1,t) ; coeff_mat_b_AL(1,t)]+...
                   [coeff_mat_k_AL(2,t) coeff_mat_k_AL(3,t);coeff_mat_b_AL(2,t) coeff_mat_b_AL(3,t)]* growth_q11_f_k;

               growth_q1_f(1,t) = coeff_mat_y_AL(1,t) + [coeff_mat_y_AL(2,t) coeff_mat_y_AL(3,t)]*[Dlogk(1,t); Dlogb(1,t)]; 
               growth_q2_f(1,t) = coeff_mat_y_AL(1,t) + [coeff_mat_y_AL(2,t) coeff_mat_y_AL(3,t)]*growth_q2_f_k; 
               growth_q3_f(1,t) = coeff_mat_y_AL(1,t) + [coeff_mat_y_AL(2,t) coeff_mat_y_AL(3,t)]*growth_q3_f_k;
               growth_q4_f(1,t) = coeff_mat_y_AL(1,t) + [coeff_mat_y_AL(2,t) coeff_mat_y_AL(3,t)]*growth_q4_f_k;
               growth_q5_f(1,t) = coeff_mat_y_AL(1,t) + [coeff_mat_y_AL(2,t) coeff_mat_y_AL(3,t)]*growth_q5_f_k;
               growth_q6_f(1,t) = coeff_mat_y_AL(1,t) + [coeff_mat_y_AL(2,t) coeff_mat_y_AL(3,t)]*growth_q6_f_k;
               growth_q7_f(1,t) = coeff_mat_y_AL(1,t) + [coeff_mat_y_AL(2,t) coeff_mat_y_AL(3,t)]*growth_q7_f_k;
               growth_q8_f(1,t) = coeff_mat_y_AL(1,t) + [coeff_mat_y_AL(2,t) coeff_mat_y_AL(3,t)]*growth_q8_f_k;
               growth_q9_f(1,t) = coeff_mat_y_AL(1,t) + [coeff_mat_y_AL(2,t) coeff_mat_y_AL(3,t)]*growth_q9_f_k;
               growth_q10_f(1,t) = coeff_mat_y_AL(1,t) + [coeff_mat_y_AL(2,t) coeff_mat_y_AL(3,t)]*growth_q10_f_k;
               growth_q11_f(1,t) = coeff_mat_y_AL(1,t) + [coeff_mat_y_AL(2,t) coeff_mat_y_AL(3,t)]*growth_q11_f_k;
               growth_q12_f(1,t) = coeff_mat_y_AL(1,t) + [coeff_mat_y_AL(2,t) coeff_mat_y_AL(3,t)]*growth_q12_f_k;
       end
       
   
    end % end Tend loop over simul periods
    growth_fore =  growth_q1_f(8)+growth_q2_f(8)+ growth_q3_f(8)+growth_q4_f(8)...
        +growth_q5_f(8)+growth_q6_f(8)+growth_q7_f(8)+growth_q8_f(8)+...
        +growth_q9_f(8)+growth_q10_f(8)+growth_q11_f(8)+growth_q12_f(8); 
    cum_grow = sum(DlogY(9:20));
    growth_fore_err =  cum_grow- growth_fore;
    
    
        
    Yg2008 = sum(DlogY(1:4));
    Yg2009 = sum(DlogY(5:8));
    Yg2010 = sum(DlogY(9:12));
    Yg2011 = sum(DlogY(13:16));
    Yg2012 = sum(DlogY(17:20));
    mean_growth_0712 = mean([4*log(gamma) Yg2008 Yg2009 Yg2010 Yg2011 Yg2012]);
    mean_growth_1316 = 4*mean(DlogY(21:36));
    
    if kkk == 1
    display('change in potential growth from 2012Q4 to 2016Q3')
    4*(y_tg(35)-y_tg(20))

    display('average GDP growth 2013 -- 2016 minus 2007 -- 2012 ')
    mean_growth_1316- mean_growth_0712
    
    display('severity of recession')
    
    mean([Yg2008 Yg2009 Yg2010])-log(gamma)*4
    
     display('average revision: (mean(y_tg(12)) - log(gamma))*4')
     (mean(y_tg(12)) - log(gamma))*4

    end

     [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

           %%% 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; 

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 = 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
              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; 
            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;
            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;
            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]; 
               [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
           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 j==1
              bbRE(j) =   -by_ratio* (1+btRE(j)-1/(gamma*exp(shockmat(j))))/(1+yyRE(j));
              
          else
                     bbRE(j) =   -by_ratio* (1+btRE(j)-(1+btRE(j-1))/(gamma*exp(shockmat(j))))/(1+yyRE(j));
          end
   end
   


PrRE =  yyRE  - HRE; 


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; 
  
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 =semi_ela*xtg*100; 
elseif kkk == 2
    trg_without = y_tg*400; % without austerity
    og_bias_without =-xtg*100; 
    cab_bias_without =semi_ela*xtg*100; 
end

if kkk == 1
display('change in CAB: 2009Q4 -- 2012Q4')
cab_with(20)-cab_with(8)

display('change in potential output growth: 2009Q4 -- 2012Q4')
trg_with(20)-trg_with(8)

display('growth_fore_err')
growth_fore_err

end
 end