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

clear all;
close all;
clc;
delete mydiary1.txt 
diary mydiary1.txt
diary on
Font_size = 11;  
Line_width = 1.5; 
load coeff_k;
load coeff_u; 
display('start')

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

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


delta = 0.084/4;  %% depreciation rate

alpha = 0.39;   %0.34; %% capital share

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

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

%tauk = 0.36;
tauk = 0.42;

%tauH = 0.4;
tauH = 0.49;

tauC = 0.17;  % 0.17
gy_ratio = 0.24; 
 
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 =  1.07;%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;

lamda3 = 1.1;

semi_ela = 0.55; 

%sigma_x = 0.0048;
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_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_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 = 100; % no. of repetitions

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

% gain parameters 

gain_control =1  ;% 1: constant gain 0 : decreasing gain
if gain_control ==1 
    

 gaink = 0.003; 
 gainb = 0.003; 
 gainw = 0.003;% 0.011;%0.004; %0
 gainR = 0.015;% 0.011;%0.005; %0
 gainkm = 0.003; 
 gainy = 0.003; %  does not matter for the temporary equilibrium
 gainy1 = gainw +0.032; 
%  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

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

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

%  xx = [0 0 0 0 0];
%  main_calibration_km_IRF(xx);
 
 xx = [gainw gainR gaink gainkm gainb];  
  if shockcontrol == 1   %   IRF following one-off shock
      main_calibration_km_IRF(xx);
  elseif shockcontrol == 2 % IRF following a sequence of shock and fiscal consolidation
      main_calibration_km_IRF_great_recession(xx);
  elseif  shockcontrol == 3                  % IRF following a sequence of shocks and fiscal consolidation
       main_calibration_km_IRF_great_recession_fc(xx);
  end
 else 
  
   
 end
display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++');
display('RE solution')
display('ordering: rk wage output mp hours cons invst caput cap gama')
display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++');
diary off

[coeff_vec_k_RE coeff_vec_s_RE]


%

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);
coeff_mat_y_AL = zeros(3, nsim);


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

% Record DlogR, DlogW, the dependent variables

Dlogw = zeros(n1, nsim);   

Dlogw0 = log(gamma);

DlogR = zeros(n1, nsim);

DlogR0 = 0; 

Dlogk = zeros(n1, nsim);

Dlogk0 = log(gamma);
Dlogkm1 = log(gamma); 

Dlogb = zeros(n1, nsim);

Dlogb0 = log(gamma);
Dlogbm1 = log(gamma); 

Dlogy0 = log(gamma); 

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

x_vec_m1 = [1;Dlogkm1; Dlogbm1]; 
x_vec_ini = [1; Dlogk0; Dlogb0]; 
%%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_tg0 = log(gamma);
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); 

% intialize forecast errors (FE)
% ForecastErrW = zeros(n1, nsim);   % FE in W equation??
 ForecastErrR = zeros(n1, nsim);   % FE in R equation??
% ForecastErrY = zeros(n1, nsim);   % FE in R equation??
% 
 ForecastErrRlevel1Q = zeros(n1, nsim);   % ???
% ForecastErrYlevel1Q = zeros(n1, nsim);

% pf_counter = zeros(n1,nsim); 
% 
% proj_u = 1.0;
% proj_l = 0.985; 
% 
% Trendgrowthest = zeros(n1, nsim); 
% Trendgrowthest_pr = zeros(n1, nsim); 
% %options = optimoptions('fsolve','Display','off');  % Turn off display
% coeff_mat_k_AL_inter = zeros(n1, nsim);
% coeff_mat_k_AL_slop = zeros(n1,nsim); 
% 
% coeff_mat_R_AL_inter = zeros(n1, nsim);
% coeff_mat_R_AL_slop = zeros(n1,nsim); 
% 
% coeff_mat_w_AL_inter = zeros(n1, nsim);
% coeff_mat_w_AL_slop = zeros(n1,nsim); 

% coeff_mat_k_AL1999 = zeros(n1,3);
% coeff_mat_b_AL1999 = zeros(n1,3);
% coeff_mat_R_AL1999 = zeros(n1,3); 
% coeff_mat_w_AL1999 = zeros(n1,3); 
% 
% y_tg1999 = zeros(n1,1); 

coeff_mat_k_AL2000 = zeros(n1,3);
coeff_mat_b_AL2000 = zeros(n1,3);
coeff_mat_R_AL2000 = zeros(n1,3); 
coeff_mat_w_AL2000 = zeros(n1,3); 


y_tg2000 = zeros(n1,1); 

Dlogy2000 = zeros(n1,1);
Dlogw2000 = zeros(n1,1);
DlogR2000 = zeros(n1,1);

Dlogk1999 = zeros(n1,1);
Dlogb1999 = zeros(n1,1);

Dlogk2000 = zeros(n1,1);
Dlogb2000 = zeros(n1,1);
khat2000 = zeros(n1,1);
bhat2000 = zeros(n1,1);
what2000 = zeros(n1,1);
RKhat2000 = zeros(n1,1);
yhat2000 = zeros(n1,1);
cab1996 =zeros(n1,1);
cab1997 =zeros(n1,1);
cab1998 =zeros(n1,1);
cab1999 =zeros(n1,1);
cab2000 =zeros(n1,1);

pg_revision = zeros(n1,nsim); % revision in potential growth rates
% growth_fe = zeros(n1, nsim); % forecast errors in output growth
% growth_fe4 = zeros(n1, nsim); % 1-year ahead forecast errors in output growth

growth_q1 = zeros(n1, nsim); % growth rate
% growth_q2 = zeros(n1, nsim); % growth rate in Q2
% growth_q3 = zeros(n1, nsim); % growth rate in Q3
% growth_q4 = zeros(n1, nsim); % growth rate in Q4

growth_annual = zeros(n1, 11); 


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); % 4Q-ahead growth forecast
growth_q6_f = zeros(n1, nsim); % 4Q-ahead growth forecast
growth_q7_f = zeros(n1, nsim); % 4Q-ahead growth forecast
growth_q8_f = zeros(n1, nsim); % 4Q-ahead growth forecast
growth_fe= zeros(n1, nsim); % forecast errors

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_annual_f =zeros(n1,11); 


AA = zeros(2,2);
BB = zeros(2,1);
mukb = zeros(2,1); 

og_shat = zeros(n1,nsim); %subjective output gaps
cab = zeros(n1,nsim);  % structural balances
bbhat = zeros(n1,nsim); %budget balance
gshock_vec = zeros(n1,nsim);

for i = 1:n1 % loop over repetions
    
    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_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);
           % 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) ;
           
           
           coeff_mat_y_AL(:,t) = coeff_mat_y_AL_ini + gainy*x_vec_m1*(Dlogy0-coeff_mat_y_AL_ini'*x_vec_m1) ;
%            
%            coeff_mat_pr_AL(:,t) = coeff_mat_pr_AL_ini + gainy*x_vec_m1*(Dlogpr0-coeff_mat_pr_AL_ini'*x_vec_m1) ;
           
          % ini_con = [shockmat(i,t) khat1]; 
           y_tg(i,t) = y_tg0; 
 %          growth_q1(i,t) = Dlogy0;
           % forecast error for wage growth
     %      ForecastErrW(i,t) = w0-wmin1+log(gamma)+shock0-psiw0; 
           % forecast error for interest rate growth
         %ForecastErrR(i,t) = DlogR0-coeff_mat_R_AL_ini'*x_vec_m1; 
         ForecastErrR(i,t) = DlogR0-coeff_mat_R_AL_ini'*x_vec_m1; 
     %     pg_revision(i,t) = 0; 
 %         growth_fe(i,t) = Dlogy0-coeff_mat_y_AL_ini'*x_vec_m1;
         
           
       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);
           
           coeff_mat_y_AL(:,t) = coeff_mat_y_AL(:,t-1) + gainy*x_vec_ini*(Dlogy(i,t-1)-coeff_mat_y_AL(:,t-1)'*x_vec_ini) ; 
%            
%            coeff_mat_pr_AL(:,t) = coeff_mat_pr_AL(:,t-1) + gainy*x_vec_ini*(Dlogpr(i,t-1)-coeff_mat_pr_AL(:,t-1)'*x_vec_ini) ; 
%            AA = [coeff_mat_b_AL(2,t) -(1-coeff_mat_b_AL(3,t)); (1-coeff_mat_k_AL(2,t)) -coeff_mat_k_AL(3,t)];
%            BB = [-coeff_mat_b_AL(1,t); coeff_mat_k_AL(1,t)];
%            mukb = AA\BB; 
%            y_tg(i,t) = coeff_mat_y_AL(1,t) + coeff_mat_y_AL(2:3,t)'*mukb; 
%            pg_revision(i,t) = y_tg(i,t) -y_tg(i,t-1);          
           
           y_tg(i,t) = y_tg(i,t-1)+gainy1*(Dlogy(i,t-1)-y_tg(i,t-1)); 
   %        growth_q1(i,t) = Dlogy(i,t-1);
            
   %        pg_revision(i,t) = y_tg(i,t) -y_tg(i,t-1); 
  %         growth_fe(i,t) = Dlogy(i,t-1)-coeff_mat_y_AL(:,t-1)'*x_vec_ini;
          % ini_con = [shockmat(i,t) khat(i,t-1)]; 
           % forecast error for wage growth
     %      ForecastErrW(i,t) = w0-wmin1+log(gamma)+shock0-psiw0; 
           % forecast error for interest rate growth
  %         ForecastErrR(i,t) = DlogR0-coeff_mat_R_AL_ini'*x_vec_m1; 
%            ForecastErrW(i,t) = what(i,t-1)-w0+log(gamma)+shockmat(i,t-1)...
%                -psiw(i,t-1);
         %   ForecastErrR(i,t) = RKhat(i,t-1)-r0-psir(i,t-1); 
          ForecastErrR(i,t) =  DlogR(i,t-1)-coeff_mat_R_AL(:,t-1)'*x_vec_ini; 
       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);
        %   ini_con = [shockmat(i,t) khat(i,t-1)]; 


%            if (coeff_mat_k_AL(2,t)>=proj_l)
%                coeff_mat_k_AL(:,t) = proj_l+(coeff_mat_k_AL(:,t-1)-proj_l)*(proj_u-proj_l)/(coeff_mat_k_AL(:,t-1)+proj_u-2*proj_l) ;
%                Mom_mat_k(:, 1 + 2*(t-1):2*t) = Mom_mat_k(:, 1 + 2*(t-2):2*(t-1));
%            end
%            if coeff_mat_k_AL(2,t)>=1
%                coeff_mat_k_AL(:,t) = coeff_mat_k_AL(:,t-1) ;
%                Mom_mat_k(:, 1 + 2*(t-1):2*t) = Mom_mat_k(:, 1 + 2*(t-2):2*(t-1));
%            else
%                if sumc(y_vec1>0.15)>=1 
%                   coeff_mat_k_AL(:,t) = coeff_mat_k_AL(:,t-1) ;
%                   Mom_mat_k(:, 1 + 2*(t-1):2*t) = Mom_mat_k(:, 1 + 2*(t-2):2*(t-1));                
%                end
%            end
%            
%            
%            if 
               
           % 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);
           
           coeff_mat_y_AL(:,t) = coeff_mat_y_AL(:,t-1) + gainy*x_vec*(Dlogy(i,t-1)-coeff_mat_y_AL(:,t-1)'*x_vec) ; 
           y_tg(i,t) = y_tg(i,t-1)+gainy1*(Dlogy(i,t-1)-y_tg(i,t-1)); 

%             AA = [coeff_mat_b_AL(2,t) -(1-coeff_mat_b_AL(3,t)); (1-coeff_mat_k_AL(2,t)) -coeff_mat_k_AL(3,t)];
%            BB = [-coeff_mat_b_AL(1,t); coeff_mat_k_AL(1,t)];
%            mukb = AA\BB; 
%            
%            
%    
%            y_tg(i,t) = coeff_mat_y_AL(1,t) + coeff_mat_y_AL(2:3,t)'*mukb; 
%            pg_revision(i,t) = y_tg(i,t) -y_tg(i,t-1);
 %           pg_revision(i,t) = gainy1*(Dlogy(i,t-1)-y_tg(i,t-1));

 %           growth_fe(i,t) =  Dlogy(i,t-1)-coeff_mat_y_AL(:,t-1)'*x_vec;
%            
%            coeff_mat_pr_AL(:,t) = coeff_mat_pr_AL(:,t-1) + gainy*x_vec*(Dlogpr(i,t-1)-coeff_mat_pr_AL(:,t-1)'*x_vec) ; 
           
%            if (coeff_mat_k_AL(2,t)>=0.995)|(coeff_mat_k_AL(2,t)<=0)
%                pf_counter(i,t) = 1; 
%      %          x_vec = [1; Dlogk(i,t-3)];
%                coeff_mat_k_AL(:,t) = coeff_mat_k_AL(:,t-1) ;
%                coeff_mat_R_AL(:,t) = coeff_mat_R_AL(:,t-1);
%                coeff_mat_w_AL(:,t) = coeff_mat_w_AL(:,t-1);
%                coeff_mat_y_AL(:,t) = coeff_mat_y_AL(:,t-1);
%                coeff_mat_pr_AL(:,t) = coeff_mat_pr_AL(:,t-1);
%            end                     

%            psiw(i,t) = psiw(i,t-1) +gainw*(what(i,t-1)-what(i,t-2)+log(gamma)...
%                +shockmat(i,t-1)- psiw(i,t-1)); 
      %     psir(i,t) = psir(i,t-1) + gainr*(RKhat(i,t-1)-RKhat(i,t-2)-psir(i,t-1)); 
%     %       ForecastErrW(i,t) = what(i,t-1)-what(i,t-2)+log(gamma)+shockmat(i,t-1)...
%                -psiw(i,t-1);
         %ForecastErrR(i,t) = RKhat(i,t-1)-RKhat(i,t-2)-psir(i,t-1);            
         ForecastErrR(i,t) = DlogR(i,t-1)-coeff_mat_R_AL(:,t-1)'*x_vec;
       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); 
          %xtg(i,t) = (t+16)*(y_tg(i,t)-log(gamma)); 
    xtg(i,t) = t*(y_tg(i,t)-log(gamma)); 
%            if t==1
%               ini_con = [shockmat(i,t) khat1]; 
%            else
%               ini_con = [shockmat(i,t) khat(i,t-1)]; 
%            end
           %      xw : trend belief in w equation
           %      xR : trend belief in R equation
           % ini_con: ordering gam lagk
      if t<=5
                 gshock = 0;                     
      elseif t == 6
                 if cab(i,t-5)<-0.005 
                    gshock = lamda3*(cab(i,t-5)+0.005);       
                 else 
                    gshock = 0;
                 end

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

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

      end
           gshock_vec(i,t) = gshock; 
           if t==1 
               ini_con = [shockmat(i,t) khat1 bhat1 xtg(i,t) gshock]; 
           else
               ini_con = [shockmat(i,t) khat(i,t-1) bhat(i,t-1) xtg(i,t) gshock]; 
           end
           
           [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
           %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
%             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(i,t) = yhat(i,t) - xtg(i,t);
           cab(i,t) = bbhat(i,t) - semi_ela* og_shat(i,t); % structural balance ratio
%         display('i')
%         i
%         display('t')
%         t
%        
% 
%            display('gshock')
%            gshock
%            

%            if t<=nsim-1
%                khat(i,t+1)=ik_ratio*ihat(i,t)+(1-delta)/gamma*(khat(i,t)-...
%                    shockmat(i,t));
%            end

          if t==1
              
           Dlogw(i,t) = what(i,t)-0+log(gamma)+shockmat(i,t); 
           DlogR(i,t) = RKhat(i,t)-0; 
           Dlogk(i,t) = khat(i,t)-0+log(gamma)+shockmat(i,t);   
           Dlogy(i,t) = yhat(i,t) - 0 + log(gamma)+shockmat(i,t);   
           Dlogb(i,t) = bhat(i,t) - 0 + 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
            if t==1 
                ForecastErrRlevel1Q(i,t) = -DlogR0;
                ForecastErrYlevel1Q(i,t) = -Dlogy0;
%                
            elseif t==2
                ForecastErrRlevel1Q(i,t) = DlogR(i,t-1)-coeff_mat_R_AL(:,t-1)'*x_vec_ini;
                ForecastErrYlevel1Q(i,t) = Dlogy(i,t-1)-coeff_mat_y_AL(:,t-1)'*x_vec_ini;
            else
                ForecastErrRlevel1Q(i,t) = DlogR(i,t-1)-coeff_mat_R_AL(:,t-1)'*x_vec; 
                ForecastErrYlevel1Q(i,t) = Dlogy(i,t-1)-coeff_mat_y_AL(:,t-1)'*x_vec; 
            end
%     coeff_mat_k_AL_inter(i,t) = coeff_mat_k_AL(1,t);
%     coeff_mat_k_AL_slop(i,t) = coeff_mat_k_AL(2,t);
%     
%     coeff_mat_R_AL_inter(i,t) = coeff_mat_R_AL(1,t);
%     coeff_mat_R_AL_slop(i,t) = coeff_mat_R_AL(2,t);
%     
%     coeff_mat_w_AL_inter(i,t) = coeff_mat_w_AL(1,t);
%     coeff_mat_w_AL_slop(i,t) = coeff_mat_w_AL(2,t);    
           if  t==nsim--4 
               cab1996(i) = cab(i,t);
           elseif t== nsim -3
                cab1997(i) = cab(i,t);
           elseif t== nsim-2
                cab1998(i) = cab(i,t);
           elseif t== nsim-1
                cab1999(i) = cab(i,t);
                Dlogk1999(i) = Dlogk(i,t); 
                Dlogb1999(i) = Dlogb(i,t); 
            elseif t==nsim
               cab2000(i) = cab(i,t);
               coeff_mat_k_AL2000(i,:) = coeff_mat_k_AL(:,t)';
               coeff_mat_b_AL2000(i,:) = coeff_mat_b_AL(:,t)';
               coeff_mat_R_AL2000(i,:) = coeff_mat_R_AL(:,t)';
               coeff_mat_w_AL2000(i,:) = coeff_mat_w_AL(:,t)';
               y_tg2000(i) = y_tg(i,t); 
               Dlogy2000(i) = Dlogy(i,t); 
               Dlogw2000(i) = Dlogw(i,t); 
               DlogR2000(i) = DlogR(i,t); 
               Dlogk2000(i) = Dlogk(i,t); 
               Dlogb2000(i) = Dlogb(i,t);                
               khat2000(i) = khat(i,t);
               bhat2000(i) = bhat(i,t);
               what2000(i) = what(i,t);
               RKhat2000(i) = RKhat(i);
               yhat2000(i) = yhat(i);

               

           end    
               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(i,t); Dlogb(i,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_q1_f(i,t) = coeff_mat_y_AL(1,t) + [coeff_mat_y_AL(2,t) coeff_mat_y_AL(3,t)]*[Dlogk(i,t); Dlogb(i,t)]; 
               growth_q2_f(i,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(i,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(i,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(i,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(i,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(i,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(i,t) = coeff_mat_y_AL(1,t) + [coeff_mat_y_AL(2,t) coeff_mat_y_AL(3,t)]*growth_q8_f_k;
               
               growth_q1(i,t) = Dlogy(i,t);
               if t >7
                  growth_fe(i,t) = growth_q5_f(i,t-4)+growth_q6_f(i,t-5)+growth_q7_f(i,t-6)+growth_q8_f(i,t-7)-growth_q1(i,t)-growth_q1(i,t-1)-growth_q1(i,t-2)-growth_q1(i,t-3); 
               end     
    end % end Tend loop over simul periods

end % end n1 loop over repetition

autocorr_FEYlevel4Q  = zeros(n1,1);  % 1st order autocorrelation for real gdp growth 4Q forecast errors
for i = 1:n1
    autocorr_FEYlevel4Q(i) = corr(growth_fe(i,1:end-1)', growth_fe(i,2:end)'); 
end
ave_autocorr_FEYlevel4Q = mean(autocorr_FEYlevel4Q);


% display('correlation between growth forecast errors and revision in potential growth rates')
% corr_i = zeros(n1,1);
% for i = 1:n1
%    corr_i(i)= corr(pg_revision(i,2000:end)',growth_fe(i,2000:end)');
% end
% median(corr_i)
% % Initial conditions
% 
% display('correlation between growth forecast errors and revision in potential growth rates')
% corr_i1 = zeros(n1,1);
% for i = 1:n1
%    corr_i1(i)= corr(pg_revision(i,2000:end-1)',growth_fe(i,2001:end)');
% end
% median(corr_i1)

display('correlation between growth forecast errors and revision in potential growth rates')
corr_i1 = zeros(n1,1);

% revised to use semi-annual frequency
growth_fe1=zeros(n1, 12); % spring to autumn vintage 


pg_revision1 = zeros(n1,12); % spring to autumn vintage 

growth_fe2=zeros(n1, 12); % autumn to spring vintage 


pg_revision2 = zeros(n1,12); % autumn to spring vintage 

% @ start of deactivation

for i =1:n1   
    for j =1:12
        % Spring to autumn vintage 
%        if mod(j,2)==0  % j is even;  Autumn vintage
% display('i')
% i
% display('j')
% j
      if j==1
               growth_fe1(i,j) = sum(growth_q1(i,(j-1)*4+1:(j-1)*4+2)) + growth_q1_f(i,(j-1)*4+2) + growth_q2_f(i,(j-1)*4+2) -   ...
               4*log(gamma);
           pg_revision1(i,j) = y_tg(i,j*4-2) - log(gamma); 
        % Autumn to spring vintage 
%        else  % j is odd; Spring vintage 
            growth_fe2(i,j) = sum(growth_q1(i,(j-1)*4+3:(j-1)*4+4)) -   ...
                growth_q1_f(i,(j-1)*4+2)    - growth_q2_f(i,(j-1)*4+2);
            pg_revision2(i,j) = y_tg(i,(j-1)*4+4) - y_tg(i,(j-1)*4+2); 
      else    
           growth_fe1(i,j) = sum(growth_q1(i,(j-1)*4+1:(j-1)*4+2)) + growth_q1_f(i,(j-1)*4+2) + growth_q2_f(i,(j-1)*4+2) -   ...
               growth_q1_f(i,(j-1)*4)    - growth_q2_f(i,(j-1)*4) - growth_q3_f(i,(j-1)*4) - growth_q4_f(i,(j-1)*4);
           pg_revision1(i,j) = y_tg(i,j*4-2) - y_tg(i,j*4-4); 
        % Autumn to spring vintage 
%        else  % j is odd; Spring vintage 
            growth_fe2(i,j) = sum(growth_q1(i,(j-1)*4+3:(j-1)*4+4)) -   ...
                growth_q1_f(i,(j-1)*4+2)    - growth_q2_f(i,(j-1)*4+2);
            pg_revision2(i,j) = y_tg(i,(j-1)*4+4) - y_tg(i,(j-1)*4+2); 
      end
        %end 
    end
end

% @ end of deactivation
% pg_revision_annual = zerons(n1, 11); % 11 year data
% growth_fe = zerons(n1, 11); % 11 year data
% for i = 1:n1
%    for j =1:11 
%       pg_revision_annual(i,j)   = pg_revision(i,2000+4*(j-1));
%    end
% 
% end

% @ start of deactivation
growth_fe = [growth_fe1 growth_fe2]; 

pg_revision = [pg_revision1 pg_revision2];

corr_growth_fe1=zeros(n1,1);
corr_growth_fe = zeros(n1,1);
for i = 1:n1
   corr_i1(i)= corr(pg_revision(i,:)',growth_fe(i,:)');
   corr_growth_fe(i) = corr(growth_fe(i,1:end-1)',growth_fe(i,2:end)');
end
median(corr_growth_fe)
std(corr_growth_fe)


% @ end of deactivation

% display('correlation between growth forecast errors and revision in potential growth rates')
% corr_i1 = zeros(n1,1);
% 
% growth_fe=zeros(n1, 11);
% pg_revision = zeros(n1,11);
% for i =1:n1
%     for j =1:11
%         growth_fe(i,j) = sum(growth_q1(i,2000+(j-1)*4+1:2000+j*4)) -   ...
%             growth_q1_f(i,2000+(j-1)*4)    -growth_q2_f(i,2000+(j-1)*4)-growth_q3_f(i,2000+(j-1)*4)-growth_q4_f(i,2000+(j-1)*4);
%         pg_revision(i,j) = y_tg(i,2000+j*4) - y_tg(i,2000+(j-1)*4); 
%     end
% end
% % pg_revision_annual = zerons(n1, 11); % 11 year data
% % growth_fe = zerons(n1, 11); % 11 year data
% % for i = 1:n1
% %    for j =1:11 
% %       pg_revision_annual(i,j)   = pg_revision(i,2000+4*(j-1));
% %    end
% % 
% % end
% for i = 1:n1
%    corr_i1(i)= corr(pg_revision(i,:)',growth_fe(i,:)');
% end
% median(corr_i1)



HRE = zeros(n1, Tend);
yy  = zeros(n1, Tend);
cRE  = zeros(n1, Tend);
wRE = zeros(n1, Tend);
REinv = zeros(n1, Tend);
PrRE = zeros(n1, Tend);
kRE = zeros(n1,Tend);
RkRE = zeros(n1,Tend);
gtRE = zeros(n1,Tend);
btRE =  zeros(n1,Tend);

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

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:n1
   for i= 1:nsim
       if i==1    
            HRE(j,i) = coeff_vec_s_RE(hours)*shockmat(j,i);
            yy(j,i) =coeff_vec_s_RE(output)*shockmat(j,i) ; 
            cRE(j,i) =coeff_vec_s_RE(cons)*shockmat(j,i);
            caputRE(j,i) = coeff_vec_s_RE(caput)*shockmat(j,i) ; 
            wRE(j,i) =coeff_vec_s_RE(wage)*shockmat(j,i);
            REinv(j,i) = coeff_vec_s_RE(invst)*shockmat(j,i);
            kRE(j,i) = coeff_vec_s_RE(cap)*shockmat(j,i);
            gtRE(j,i) = coeff_vec_s_RE(gt)*shockmat(j,i);
            btRE(j,i) = coeff_vec_s_RE(bt)*shockmat(j,i);
            RkRE(j,i) = coeff_vec_s_RE(rk)*shockmat(j,i);
            ini_con_RE = [shockmat(j,i) 0 0];
            [y_vec1_RE, fore_dis_RE] = Learningsol(xw_RE, xR_RE, xk_RE, xb_RE, ini_con_RE);
            SwRE(j,i) = fore_dis_RE(1);
            SrRE(j,i) = fore_dis_RE(2);
          
            
       else
            HRE(j,i) = coeff_vec_k_RE(hours,1)*kRE(j,i-1) + coeff_vec_k_RE(hours,2)*btRE(j,i-1) + coeff_vec_s_RE(hours)*shockmat(j,i);
            yy(j,i) = coeff_vec_k_RE(output,1)*kRE(j,i-1) + coeff_vec_k_RE(output,2)*btRE(j,i-1) + coeff_vec_s_RE(output)*shockmat(j,i); 
            cRE(j,i) =coeff_vec_k_RE(cons,1)*kRE(j,i-1) +coeff_vec_k_RE(cons,2)*btRE(j,i-1)+ coeff_vec_s_RE(cons)*shockmat(j,i);
            kRE(j,i) = coeff_vec_k_RE(cap,1)*kRE(j,i-1)+coeff_vec_k_RE(cap,2)*btRE(j,i-1) + coeff_vec_s_RE(cap)*shockmat(j,i);
            
            wRE(j,i) =coeff_vec_k_RE(wage,1)*kRE(j,i-1)+coeff_vec_k_RE(wage,2)*btRE(j,i-1) +coeff_vec_s_RE(wage)*shockmat(j,i);
            REinv(j,i) = coeff_vec_k_RE(invst,1)*kRE(j,i-1)+coeff_vec_k_RE(invst,2)*btRE(j,i-1) + coeff_vec_s_RE(invst)*shockmat(j,i);
            RkRE(j,i) = coeff_vec_k_RE(rk,1)*kRE(j,i-1) +coeff_vec_k_RE(rk,2)*btRE(j,i-1)+ coeff_vec_s_RE(rk)*shockmat(j,i);
            gtRE(j,i) = coeff_vec_k_RE(gt,1)*kRE(j,i-1) +coeff_vec_k_RE(gt,2)*btRE(j,i-1)+ coeff_vec_s_RE(gt)*shockmat(j,i);    
            btRE(j,i) = coeff_vec_k_RE(bt,1)*kRE(j,i-1) +coeff_vec_k_RE(bt,2)*btRE(j,i-1)+ coeff_vec_s_RE(bt)*shockmat(j,i);
            ini_con_RE = [shockmat(j,i) kRE(j,i-1) btRE(j,i-1)];
            [y_vec1_RE, fore_dis_RE] = Learningsol(xw_RE, xR_RE, xk_RE, xb_RE, ini_con_RE);
            SwRE(j,i) = fore_dis_RE(1);
            SrRE(j,i) = fore_dis_RE(2);            
           
       end

   end
end

PrRE=  yy - HRE; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Plot the RE & Learning dynamics
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%kk = [zeros(n1,1) kk(:, 1:end-1)];    % add the initial predetrmined period
figure 

% only median of the simulation is plotted
subplot(3,3,1)
set(gca,'fontsize',Font_size)
% capital stock plot
plot(median(kRE,1), 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(median(khat,1),'r-+')
title('Capital')

%output plot
subplot(3,3,2)
set(gca,'fontsize',Font_size)
plot(median(yy,1), 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(median(yhat,1),'r-+')
title('Output')

% investment plot
subplot(3,3,3)
set(gca,'fontsize',Font_size)
plot(median(REinv,1), 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(median(ihat,1),'r-+')
title('Investment')

% consumption plot
subplot(3,3,4)
set(gca,'fontsize',Font_size)
plot(median(cRE,1), 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(median(chat,1),'r-+')
title('Consumption')

% hours plot
subplot(3,3,5)
set(gca,'fontsize',Font_size)
plot(median(HRE,1), 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(median(Hhat,1),'r-+')
title('Hours')

% wage plot
subplot(3,3,6)
set(gca,'fontsize',Font_size)
plot(median(wRE,1), 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(median(what,1),'r-+')
title('Wage')

% rental rate plot
subplot(3,3,7)
set(gca,'fontsize',Font_size)
plot(median(RkRE,1), 'b-o')
hold on 
plot(zeros(Tend,1))
hold on 
plot(median(RKhat,1),'r-+')
title('Rental Rate')

% income sum plot
subplot(3,3,8)
set(gca,'fontsize',Font_size)
%plot(median(SyRE,1), 'b-o')
%hold on 
plot(zeros(Tend,1))
hold on 
plot(median(Swe,1),'r-+')
title('Infinite Sum of Wage(AL)/Income')

% interest rate sum plot
subplot(3,3,9)
set(gca,'fontsize',Font_size)
%plot(median(SrRE,1), 'b-o')
%hold on 
plot(zeros(Tend,1))
hold on 
plot(median(Sre,1),'r-+')
title('Infinite Sum of Interest Rate')
           
% figure 
% % evolution of beliefs : for all simulations
% subplot(2,1,1)
% set(gca,'fontsize',Font_size)
% plot(psiw')
% hold on
% plot(log(gamma)*ones(1,Tend))
% title('Evolution of \psi w')
% subplot(2,1,2)
% plot(psir')
% hold on
% plot(zeros(1,Tend))
% title('Evolution of \psi r')
% 
% save psiw.mat psiw


ln_H0 = log(0.3); % set initial ss hours at 0.3
ln_Y0 = 0;     % normalization X0 such that y_bar/X0 = 1  
ln_C0 = log(cy_ratio); 
%ln_C0 = 0.1

ln_W0 = log(1-alpha)-ln_H0;  
ln_I0 = log(1-exp(ln_C0));
ln_Pr = log(1/0.3); 

ini_vec = [ln_Y0*ones(n1,1); ln_C0*ones(n1,1);  ln_W0*ones(n1,1); ln_I0*ones(n1, 1); ln_Pr*ones(n1,1)];

% by default the model starts from the steady state 

shocks_mat = [repmat(shockmat(:,1:nsim), 5,1)] ;  % note hours are stationary variables


% order H y c w I 
% Learning artificial data 

% difference of log deviations: Learning

H_diff = [Hhat(:,1) (Hhat(:,2:end)-Hhat(:,1:end-1))];  
y_diff = [yhat(:,1) (yhat(:,2:end)-yhat(:,1:end-1))];
c_diff = [chat(:,1) (chat(:,2:end)-chat(:,1:end-1))];
w_diff = [what(:,1) (what(:,2:end)-what(:,1:end-1))];
I_diff = [ihat(:,1) (ihat(:,2:end)-ihat(:,1:end-1))];
Pr_diff = [Prhat(:,1) (Prhat(:,2:end)-Prhat(:,1:end-1))];

diff_mat = [y_diff;c_diff; w_diff; I_diff; Pr_diff];   % H not included because stationary
level_mat = zeros(5*n1,Tend);     % to record the level of the 5 variables
H_level = zeros(n1,Tend);   % to record the level of hours
org_growth =zeros(5*n1,Tend);

% RE data

% difference of log deviations: RE 
% no. of repetitions * no. of simul periods
% orders H y c w inv Pr
HRE_growth = [HRE(:,1) (HRE(:,2:end)-HRE(:,1:end-1))] ;   % change of log deviations different from org_growth_RE  
yRE_growth = [yy(:,1) (yy(:,2:end)-yy(:,1:end-1))]; 
cRE_growth = [cRE(:,1) (cRE(:,2:end)-cRE(:,1:end-1))] ;
wRE_growth = [wRE(:,1) (wRE(:,2:end)-wRE(:,1:end-1))] ;
REinv_growth = [REinv(:,1) (REinv(:,2:end)-REinv(:,1:end-1))] ;
PrRE_growth = [PrRE(:,1) (PrRE(:,2:end)-PrRE(:,1:end-1))] ;

growth_RE_mat = [yRE_growth; cRE_growth; wRE_growth; REinv_growth; PrRE_growth];   % H not included because stationary

level_RE_mat = zeros(5*n1,Tend); % to record the level of the 4 variables
HRE_level = zeros(n1,Tend);   % to record level of hours

org_growth_RE = zeros(5*n1,Tend);  % Growth rates based on level of data
% y_level = zeros(n1,Tend)
% c_level = zeros(n1, Tend) 
% H_level = zeros(n1, Tend)
% w_level = zeros(n1, Tend)
% I_level = zeros(n1, Tend)

for i = 1:Tend
    if i ==1 
        % Growth rates data  
        org_growth(:,i) = diff_mat(:,i) + shocks_mat(:,i)+log(gamma)*ones(5*n1,1); 
        % levels
        level_mat(:,i) = ini_vec+ org_growth(:,i);       % level data
        H_level(:,i) = ln_H0*ones(n1,1)+H_diff(:,i);   % hours level data
        
        org_growth_RE(:,i)= growth_RE_mat(:,i) + shocks_mat(:,i)+log(gamma)*ones(5*n1,1);
        level_RE_mat(:,i) = ini_vec+ org_growth_RE(:,i);          
        HRE_level(:,i) = ln_H0*ones(n1,1)+HRE_growth(:,i);        
    else
        org_growth(:,i) = diff_mat(:,i) + shocks_mat(:,i)+log(gamma)*ones(5*n1,1);
        level_mat(:,i) = level_mat(:,i-1)+org_growth(:,i);
        H_level(:,i) = H_level(:,i-1)+H_diff(:,i);    
        
        org_growth_RE(:,i) =  growth_RE_mat(:,i) + shocks_mat(:,i)+log(gamma)*ones(5*n1,1);
        level_RE_mat(:,i) = level_RE_mat(:,i-1)+ org_growth_RE(:,i);
        HRE_level(:,i) = HRE_level(:,i-1)+HRE_growth(:,i);         
    end
end



% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Autocorrelations: forecast errors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


autocorr_FEY = zeros(n1,1);  % 1st order autocorrelation for income growth 1Q forecast errors
autocorr_FER =zeros(n1,1);  % 1st order autocorrelation for interst growth 1Q FE
autocorr_FERlevel1Q = zeros(n1,1);  % 1st order autocorrelation for interest level 1Q FE
autocorr_FEYlevel1Q = zeros(n1,1);  % 1st order autocorrelation for output level 1Q FE
% 
autocorr_FERlevel4Q = zeros(n1,1);  % 1st order autocorrelation for a year ahead interest level FE
% 
% autocorr_trendgrowthfe =zeros(n1,1);
for i = 1:n1
%   autocorr_FEW(i)=corr(ForecastErrW(i,Tend1:Tend-1)', ForecastErrW(i,Tend1+1:Tend)'); 
   autocorr_FER(i)=corr(ForecastErrR(i,Tend1:Tend-1)', ForecastErrR(i,Tend1+1:Tend)'); 
   autocorr_FERlevel1Q(i) = corr(ForecastErrRlevel1Q(i,Tend1:Tend-1)', ForecastErrRlevel1Q(i,Tend1+1:Tend)');
  autocorr_FEYlevel1Q(i) = corr(ForecastErrYlevel1Q(i,Tend1:Tend-1)', ForecastErrYlevel1Q(i,Tend1+1:Tend)');
%   autocorr_FERlevel4Q(i) = corr(ForecastErrRlevel4Q(i,Tend1:Tend-5)', ForecastErrRlevel4Q(i,Tend1+1:Tend-4)');
%    
% %   autocorr_trendgrowthfe(i)= corr(trendgrowthfe(i,Tend1:Tend-1)', trendgrowthfe(i,Tend1+1:Tend)'); 
 end
% ave_autocorr_FEW = median(autocorr_FEW);
 ave_autocorr_FER = median(autocorr_FER);
% ave_autocorr_FERlevel1Q = median(autocorr_FERlevel1Q); 
% 
 ave_autocorr_FERlevel1Q = median(autocorr_FERlevel1Q);

%ave_autocorr_FERlevel4Q = median(autocorr_FERlevel4Q); 

%ave_trendgrowthfe = median(autocorr_trendgrowthfe);

% display('++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
% display('Forecast error persistence: Trend growth ')
% ave_trendgrowthfe

%[ave_autocorr_FEY ave_autocorr_FER ] 
% display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
% display('Median forecast error persistence: GDP growth1Q, R1Q')
% display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
% %[ave_autocorr_FEW ave_autocorr_FERlevel1Q ave_autocorr_FERlevel4Q] 
% [ave_autocorr_FEYlevel1Q ave_autocorr_FERlevel1Q]
% 
% std_autocorr_FEY = std(autocorr_FEYlevel1Q);
 std_autocorr_FER = std(autocorr_FERlevel1Q);

% display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
% display('Standard deviation of forecast error persistence: GDP growth1Q, R1Q')
% display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
% [std_autocorr_FEY std_autocorr_FER] 

% only use lastT period to calculate the statistics 
hp_data_sample = [H_level; level_mat];   % levels for all variables under learning: including hours
hp_data_sample = hp_data_sample(:, Tend1:Tend); 
hp_data_sample = hp_data_sample'; 

hp_data_sample_RE = [HRE_level; level_RE_mat]; 
hp_data_sample_RE = hp_data_sample_RE(:, Tend1:Tend); 
hp_data_sample_RE = hp_data_sample_RE'; 

save hp_data_sample.mat
display('size of hp_data_sample')
size(hp_data_sample)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Filtering the data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[y_sample_trend,desvabs] = hpfilter(hp_data_sample,1600);
% number of observations * number of series

[y_sample_trend_RE,desvabs_RE] = hpfilter(hp_data_sample_RE,1600);



% figure 
% plot(y_sample_trend)


var_level_detr = hp_data_sample'-y_sample_trend';  % get detrended data

var_level_detr_RE = hp_data_sample_RE'-y_sample_trend_RE';  % get detrended data

% calculating the correlation between detrended consumption and hours

%display('correlation between detrended consumption and hours')
org_growth_mat = [H_diff; org_growth];
org_growth_RE_mat = [HRE_growth; org_growth_RE];
 
%%Compute the correlation of the trend growth forecasts with HP filtered
%%variables

corr_tg_y = zeros(n1,1);
corr_tg_c = zeros(n1,1);
corr_tg_I = zeros(n1,1);
corr_tg_H = zeros(n1,1);

corr_tg_y_pr = zeros(n1,1);
corr_tg_c_pr = zeros(n1,1);
corr_tg_I_pr = zeros(n1,1);
corr_tg_H_pr = zeros(n1,1);


corr_tg = zeros(n1,1); 
corr_tg_pr = zeros(n1,1); 
corr_tg_rea = zeros(n1,1); 

%Trendgrowthest_rea = Trendgrowthest - shockmat; 
% for i = 1: n1
%     corr_tg_y(i) = corr(var_level_detr(n1+i,:)', Trendgrowthest(i,Tend1:Tend)'); 
%     corr_tg_c(i) = corr(var_level_detr(2*n1+i,:)',    Trendgrowthest(i,Tend1:Tend)' ); 
%     corr_tg_I(i) = corr(var_level_detr(4*n1+i,:)'  ,  Trendgrowthest(i,Tend1:Tend)' );
%     corr_tg_H(i) = corr(var_level_detr(i,:)' ,   Trendgrowthest(i,Tend1:Tend)' );
%     
%     corr_tg_y_pr(i) = corr(var_level_detr(n1+i,:)', Trendgrowthest_pr(i,Tend1:Tend)'); 
%     corr_tg_c_pr(i) = corr(var_level_detr(2*n1+i,:)',    Trendgrowthest_pr(i,Tend1:Tend)' ); 
%     corr_tg_I_pr(i) = corr(var_level_detr(4*n1+i,:)'  ,  Trendgrowthest_pr(i,Tend1:Tend)' );
%     corr_tg_H_pr(i) = corr(var_level_detr(i,:)' ,   Trendgrowthest_pr(i,Tend1:Tend)' );
%     
%     
%     corr_tg(i) = corr(Trendgrowthest(i,Tend1:Tend-1)' ,   Trendgrowthest(i,Tend1+1:Tend)' );
%    % corr_tg_rea(i) = corr(Trendgrowthest_rea(i,Tend1:Tend-1)' ,   Trendgrowthest_rea(i,Tend1+1:Tend)' );
%     corr_tg_pr(i) = corr(Trendgrowthest_pr(i,Tend1:Tend-1)' ,   Trendgrowthest_pr(i,Tend1+1:Tend)' );
% end
% corr_tg_y_med = median(corr_tg_y);
% corr_tg_y_std = std(corr_tg_y);
% corr_tg_c_med = median(corr_tg_c);
% corr_tg_c_std = std(corr_tg_c);
% corr_tg_I_med = median(corr_tg_I);
% corr_tg_I_std = std(corr_tg_I);
% corr_tg_H_med = median(corr_tg_H);
% corr_tg_H_std = std(corr_tg_H);
% 
% 
% corr_tg_y_pr_med = median(corr_tg_y_pr);
% corr_tg_y_pr_std = std(corr_tg_y_pr);
% corr_tg_c_pr_med = median(corr_tg_c_pr);
% corr_tg_c_pr_std = std(corr_tg_c_pr);
% corr_tg_I_pr_med = median(corr_tg_I_pr);
% corr_tg_I_pr_std = std(corr_tg_I_pr);
% corr_tg_H_pr_med = median(corr_tg_H_pr);
% corr_tg_H_pr_std = std(corr_tg_H_pr);
% 
% corr_tg_med = median(corr_tg);
% corr_tg_std = std(corr_tg);
% corr_tg_pr_med = median(corr_tg_pr);
% corr_tg_pr_std = std(corr_tg_pr);
% 
% display('median correlation between trend growth and detrended: y, c, I ,H')

% corr_tg_y_med
% corr_tg_c_med
% corr_tg_I_med 
% corr_tg_H_med
% 
% display('median correlation between trend output per hour growth and detrended: y, c, I ,H')
% 
% corr_tg_y_pr_med
% corr_tg_c_pr_med
% corr_tg_I_pr_med 
% corr_tg_H_pr_med
% 
% 
% display('std of correlation between trend growth and detrended: y, c, I ,H')
% corr_tg_y_std
% corr_tg_c_std
% corr_tg_I_std 
% corr_tg_H_std
% 
% display('std of correlation between trend growth and detrended: y, c, I ,H')
% corr_tg_y_pr_std
% corr_tg_c_pr_std
% corr_tg_I_pr_std 
% corr_tg_H_pr_std
% 
% display('autocorrelation: trend growth forecasts')
% 
% corr_tg_med
% corr_tg_std
% 
% display('autocorrelation: trend productivity growth forecasts')
% 
% corr_tg_pr_med
% corr_tg_pr_std

%% Compute bc statistics

  % Levels 
  
  % Volatilities HP-filtered data
  
  
  STD_var_level_detr = std(var_level_detr,0,2);   % Learning data s.d.
     
  STD_var_level_detr_RE = std(var_level_detr_RE,0,2);  % RE data s.d.

  ave_std_y = median(std_A/STD_var_level_detr(n1+1:2*n1)) ;  % what is std_A?
  std_std_y = std(std_A/STD_var_level_detr(n1+1:2*n1)); 
% orders H y c w inv Pr    
ave_std_cyratio = median(STD_var_level_detr(2*n1+1:3*n1)./STD_var_level_detr(n1+1:2*n1)); 
ave_std_invyratio = median(STD_var_level_detr(4*n1+1:5*n1)./STD_var_level_detr(n1+1:2*n1)); 
ave_std_Hyratio = median(STD_var_level_detr(1:n1)./STD_var_level_detr(n1+1:2*n1)); 
ave_std_Pryratio = median(STD_var_level_detr(5*n1+1:6*n1)./STD_var_level_detr(n1+1:2*n1)); 

ave_std_cyratio_RE = median(STD_var_level_detr_RE(2*n1+1:3*n1)./STD_var_level_detr_RE(n1+1:2*n1)); 
ave_std_invyratio_RE = median(STD_var_level_detr_RE(4*n1+1:5*n1)./STD_var_level_detr_RE(n1+1:2*n1)); 
ave_std_Hyratio_RE = median(STD_var_level_detr_RE(1:n1)./STD_var_level_detr_RE(n1+1:2*n1)); 
ave_std_Pryratio_RE = median(STD_var_level_detr_RE(5*n1+1:6*n1)./STD_var_level_detr_RE(n1+1:2*n1)); 

std_std_cyratio = std(STD_var_level_detr(2*n1+1:3*n1)./STD_var_level_detr(n1+1:2*n1)); 
std_std_invyratio = std(STD_var_level_detr(4*n1+1:5*n1)./STD_var_level_detr(n1+1:2*n1)); 
std_std_Hyratio = std(STD_var_level_detr(1:n1)./STD_var_level_detr(n1+1:2*n1)); 
std_std_Pryratio = std(STD_var_level_detr(5*n1+1:6*n1)./STD_var_level_detr(n1+1:2*n1)); 

% 
display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
display('Median standard deviation of HP filtered output ')


 % s.d. ratio c, inv, H , relative to output 
 display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
 display('Median Relative Standard Deviation: c, inv, H , Pr to output (HP-detrended data) ')
 display('column 1: learning model; column 2: RE model')
 display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
 
 
  [ave_std_cyratio ave_std_cyratio_RE; 
    ave_std_invyratio  ave_std_invyratio_RE; 
   ave_std_Hyratio  ave_std_Hyratio_RE ;
    ave_std_Pryratio  ave_std_Pryratio_RE]
    
display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
display('Standard deviation of relative stds: c, inv, H , Pr to output (HP-detrended data)')
display('learning model')
display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
[std_std_cyratio std_std_invyratio std_std_Hyratio std_std_Pryratio]

 % s.d ratio: learning/RE 
 % orders H y c w inv Pr
 ave_std_H_ratio =  median(STD_var_level_detr(1:n1)./STD_var_level_detr_RE(1:n1));
 ave_std_y_ratio =  median(STD_var_level_detr(1+n1:2*n1)./STD_var_level_detr_RE(1+n1:2*n1));
 ave_std_c_ratio =  median(STD_var_level_detr(1+2*n1:3*n1)./STD_var_level_detr_RE(1+2*n1:3*n1));
 ave_std_w_ratio =  median(STD_var_level_detr(1+3*n1:4*n1)./STD_var_level_detr_RE(1+3*n1:4*n1));
 ave_std_inv_ratio =  median(STD_var_level_detr(1+4*n1:5*n1)./STD_var_level_detr_RE(1+4*n1:5*n1));
 ave_std_Pr_ratio =  median(STD_var_level_detr(1+5*n1:6*n1)./STD_var_level_detr_RE(1+5*n1:6*n1));

 display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
 display(' Median Relative standard deviation: learning/RE model;')
 display('Ordering of variables: H y c w inv Pr (HP-detrended data)')
 display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
 % orders H y c w inv Pr


[ ave_std_H_ratio; ave_std_y_ratio ; ave_std_c_ratio; 
 ave_std_w_ratio; ave_std_inv_ratio; ave_std_Pr_ratio ]% relative s.d. Learning/RE
 % correlations (HP-detrended data)
 
     %H y correlation learning
     % correlations for all repetitions: learning modal
     corr_Hy = zeros(n1,1);   
     corr_cy = zeros(n1,1);
     corr_Iy = zeros(n1,1); 
     corr_Pry = zeros(n1,1);
     corr_PrH =zeros(n1,1);
     
     % correlations for all repetitions: RE modal
     corr_Hy_RE = zeros(n1,1);
     corr_cy_RE = zeros(n1,1);
     corr_Iy_RE = zeros(n1,1); 
     corr_Pry_RE = zeros(n1,1);
     corr_PrH_RE = zeros(n1,1);
     
    % orders H y c w inv Pr   
    
     for i =1:n1
        corr_Hy(i) =  corr(var_level_detr(i,:)', var_level_detr(i+n1,:)'); 
        corr_cy(i) =  corr(var_level_detr(i+2*n1,:)', var_level_detr(i+n1,:)'); 
        corr_Iy(i) =  corr(var_level_detr(i+4*n1,:)', var_level_detr(i+n1,:)'); 
        corr_Pry(i) =  corr(var_level_detr(i+5*n1,:)', var_level_detr(i+n1,:)'); 
        corr_PrH(i) = corr(var_level_detr(i+5*n1,:)', var_level_detr(i,:)'); 
        
        corr_Hy_RE(i) =  corr(var_level_detr_RE(i,:)', var_level_detr_RE(i+n1,:)'); 
        corr_cy_RE(i) =  corr(var_level_detr_RE(i+2*n1,:)', var_level_detr_RE(i+n1,:)'); 
        corr_Iy_RE(i) =  corr(var_level_detr_RE(i+4*n1,:)', var_level_detr_RE(i+n1,:)');        
        corr_Pry_RE(i) =  corr(var_level_detr_RE(i+5*n1,:)', var_level_detr_RE(i+n1,:)');   
        corr_PrH_RE(i) =  corr(var_level_detr_RE(i+5*n1,:)', var_level_detr_RE(i,:)');  
     end
     
     ave_corr_Hy = median(corr_Hy); ave_corr_cy = median(corr_cy); 
     ave_corr_Iy = median(corr_Iy); ave_corr_Pry = median(corr_Pry); 
     ave_corr_PrH = median(corr_PrH); 
     
     
     ave_corr_Hy_RE = median(corr_Hy_RE); ave_corr_cy_RE = median(corr_cy_RE); 
     ave_corr_Iy_RE = median(corr_Iy_RE)  ;  ave_corr_Pry_RE = median(corr_Pry_RE); 
     ave_corr_PrH_RE = median(corr_PrH_RE); 
     
     display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
     display('Median correlation cy, Iy, Hy, Pry, PrH (HP-detrended data)')
     display('Column 1: learning model; Column 2: RE model')
     display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
     
     [ave_corr_cy ave_corr_cy_RE;
      ave_corr_Iy ave_corr_Iy_RE;
      ave_corr_Hy ave_corr_Hy_RE;              
     ave_corr_Pry ave_corr_Pry_RE;
     ave_corr_PrH ave_corr_PrH_RE] 
 
     std_corr_Hy = std(corr_Hy); std_corr_cy = std(corr_cy); 
     std_corr_Iy = std(corr_Iy); std_corr_Pry = std(corr_Pry); 
     std_corr_PrH = std(corr_PrH); 

      display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
     display('standard deviation of correlations cy, Iy, Hy, Pry, PrH (HP-detrended data)')
     display('learning model')
     display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')   
     
     [std_corr_cy std_corr_Iy std_corr_Hy std_corr_Pry std_corr_PrH]
     
 display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
 display('   Growth rates   statistics ');
 display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
 
 % Autocorrelations in growth rates
%  org_growth_mat = org_growth_mat(:,Tend1:Tend);
%  
% [ygrowth_trend,desvabs_grow] = hpfilter(org_growth_mat',1600);
% corr_trend_ygrowth = zeros(n1,1);
% corr_tg_tvp_trrend =zeros(n1,1);
% 
% corr_trend_ygrowth_prod = zeros(n1,1);
% corr_tg_tvp_trrend_prod =zeros(n1,1);
% for i =1 :n1
%    corr_trend_ygrowth(i) = corr(ygrowth_trend(1:end-1,n1+i), ygrowth_trend(2:end, n1+i));
%    corr_tg_tvp_trrend(i) = corr(Trendgrowthest(i,Tend1:Tend-1)'-ygrowth_trend(1:end-1,n1+i) ,   Trendgrowthest(i,Tend1+1:Tend)'- ygrowth_trend(2:end, n1+i));
%    
%    corr_trend_ygrowth_prod(i) = corr(ygrowth_trend(1:end-1,5*n1+i), ygrowth_trend(2:end, 5*n1+i));
%    corr_tg_tvp_trrend_prod(i) = corr(Trendgrowthest_pr(i,Tend1:Tend-1)'-ygrowth_trend(1:end-1,5*n1+i) ,   Trendgrowthest_pr(i,Tend1+1:Tend)'- ygrowth_trend(2:end, 5*n1+i));
%    
% end
% display('median trend output growth')
% median(corr_trend_ygrowth)   % 0.9930
% display('median trend growth forecast error')
% median(corr_tg_tvp_trrend)   % 0.9919
% display('std trend growth forecast error')
% std(corr_tg_tvp_trrend)
% 
% display('median trend output per hour growth')
% median(corr_trend_ygrowth_prod)   % 0.9930
% display('median trend output per hour growth forecast error')
% median(corr_tg_tvp_trrend_prod) 
% display('std trend output per hour growth forecast error')
% std(corr_tg_tvp_trrend_prod) 
% 
%  org_growth_RE_mat = org_growth_RE_mat(:,Tend1:Tend);
%  simulatedY= org_growth_mat(1:n1,:)'; 
% % simulatedR = RKhat(:,Tend1-1:Tend);
% 
% simulatedR = RKhat;
% % save simulatedR.mat simulatedR; 
% % save simulatedY.mat simulatedY;
% 

   %%Growth rates volatilities
   
  STD_org_growth = std(org_growth_mat,0,2); 
   
  STD_org_growth_RE = std(org_growth_RE_mat,0,2); 
  % orders H y c w inv Pr   
%   ave_H_growth_std= median(STD_org_growth(1:n1));
   ave_y_growth_std= median(STD_org_growth(n1+1:2*n1));
   
%   ave_c_growth_std= median(STD_org_growth(2*n1+1:3*n1));
%   ave_w_growth_std= median(STD_org_growth(3*n1+1:4*n1));
%   ave_inv_growth_std =median(STD_org_growth(4*n1+1:5*n1));
%   ave_Pr_growth_std =median(STD_org_growth(5*n1+1:6*n1));
  % orders H y c w inv Pr  
%   ave_H_growth_RE_std= median(STD_org_growth_RE(1:n1));
   ave_y_growth_RE_std= median(STD_org_growth_RE(n1+1:2*n1));
%   ave_c_growth_RE_std= median(STD_org_growth_RE(2*n1+1:3*n1));
%   ave_w_growth_RE_std= median(STD_org_growth_RE(3*n1+1:4*n1));
%   ave_inv_growth_RE_std =median(STD_org_growth_RE(4*n1+1:5*n1)); 
%   ave_Pr_growth_RE_std =median(STD_org_growth_RE(5*n1+1:6*n1)); 
  
  y_growth_std400 = 400*ave_y_growth_std;
  y_growth_std400_RE = 400*ave_y_growth_RE_std;

    % orders H y c w inv Pr   
  ave_std_cy_growth = median(STD_org_growth(2*n1+1:3*n1)./STD_org_growth(n1+1:2*n1)); 
   ave_std_Iy_growth = median(STD_org_growth(4*n1+1:5*n1)./STD_org_growth(n1+1:2*n1)); 
   ave_std_Hy_growth = median(STD_org_growth(1:n1)./STD_org_growth(n1+1:2*n1)); 
   ave_std_Pry_growth = median(STD_org_growth(5*n1+1:6*n1)./STD_org_growth(n1+1:2*n1)); 
   
  ave_std_cy_growth_RE = median(STD_org_growth_RE(2*n1+1:3*n1)./STD_org_growth_RE(n1+1:2*n1)); 
   ave_std_Iy_growth_RE = median(STD_org_growth_RE(4*n1+1:5*n1)./STD_org_growth_RE(n1+1:2*n1)); 
   ave_std_Hy_growth_RE = median(STD_org_growth_RE(1:n1)./STD_org_growth_RE(n1+1:2*n1)); 
   ave_std_Pry_growth_RE = median(STD_org_growth_RE(5*n1+1:6*n1)./STD_org_growth_RE(n1+1:2*n1)); 
   
  display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
  display('median relative standard deviations(growth rates data)')
  display('c I H Pr: relative to output')
  display('Column 1 learning model; column 2 RE model')
  display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
   
  [ave_std_cy_growth ave_std_cy_growth_RE;
    ave_std_Iy_growth ave_std_Iy_growth_RE;
    ave_std_Hy_growth ave_std_Hy_growth_RE;
   ave_std_Pry_growth   ave_std_Pry_growth_RE]

  std_y_growth_std= std(STD_org_growth(n1+1:2*n1));
    display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
  display('standard deviation output growth')
    display('learning model')
  display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')  
  [std_y_growth_std*400]
  
  std_std_cy_growth = std(STD_org_growth(2*n1+1:3*n1)./STD_org_growth(n1+1:2*n1)); 
   std_std_Iy_growth = std(STD_org_growth(4*n1+1:5*n1)./STD_org_growth(n1+1:2*n1)); 
   std_std_Hy_growth = std(STD_org_growth(1:n1)./STD_org_growth(n1+1:2*n1)); 
   std_std_Pry_growth = std(STD_org_growth(5*n1+1:6*n1)./STD_org_growth(n1+1:2*n1)); 
   
   display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
  display('standard deviation of relative standard deviations(growth rates data)')
  display('c I H Pr: relative to output')
  display('learning model')
  display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')  
  [std_std_cy_growth std_std_Iy_growth std_std_Hy_growth  std_std_Pry_growth]

  % orders H y c w inv Pr
ave_std_H_growth_ratio = median(STD_org_growth(1:n1)./STD_org_growth_RE(1:n1)); 
ave_std_y_growth_ratio = median(STD_org_growth(1+n1:2*n1)./STD_org_growth_RE(1+n1:2*n1)); 
ave_std_c_growth_ratio = median(STD_org_growth(1+2*n1:3*n1)./STD_org_growth_RE(1+2*n1:3*n1)); 
ave_std_w_growth_ratio = median(STD_org_growth(1+3*n1:4*n1)./STD_org_growth_RE(1+3*n1:4*n1)); 
ave_std_inv_growth_ratio = median(STD_org_growth(1+4*n1:5*n1)./STD_org_growth_RE(1+4*n1:5*n1)); 
ave_std_Pr_growth_ratio = median(STD_org_growth(1+5*n1:6*n1)./STD_org_growth_RE(1+5*n1:6*n1)); 
       

  display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
  display('median relative standard deviation ratios: learning/RE (growth rates data)')
  display('H y c w inv Pr')
  display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
  
   [ ave_std_H_growth_ratio; 
     ave_std_y_growth_ratio;
     ave_std_c_growth_ratio;
     ave_std_w_growth_ratio;
     ave_std_inv_growth_ratio;
     ave_std_Pr_growth_ratio]     % relative s.d. Learning/RE
 
% correlation 
 
     %H y correlation learning
     corr_Hy_growth = zeros(n1,1);
     corr_cy_growth = zeros(n1,1);
     corr_Iy_growth = zeros(n1,1); 
     corr_Pry_growth = zeros(n1,1);
     corr_PrH_growth = zeros(n1,1);
     
     corr_Hy_RE_growth = zeros(n1,1);
     corr_cy_RE_growth = zeros(n1,1);
     corr_Iy_RE_growth = zeros(n1,1); 
     corr_Pry_RE_growth = zeros(n1,1);
     corr_PrH_RE_growth = zeros(n1,1);
%      
     for i =1:n1
         % orders H y c w inv Pr  
        corr_Hy_growth(i) =  corr(org_growth_mat(i,:)', org_growth_mat(i+n1,:)'); 
        corr_cy_growth(i) =  corr(org_growth_mat(i+2*n1,:)', org_growth_mat(i+n1,:)'); 
        corr_Iy_growth(i) =  corr(org_growth_mat(i+4*n1,:)', org_growth_mat(i+n1,:)'); 
        corr_Pry_growth(i) =  corr(org_growth_mat(i+5*n1,:)', org_growth_mat(i+n1,:)'); 
        corr_PrH_growth(i) =  corr(org_growth_mat(i+5*n1,:)', org_growth_mat(i,:)'); 
        % orders H y c w inv Pr  
        corr_Hy_RE_growth(i) =  corr(org_growth_RE_mat(i,:)', org_growth_RE_mat(i+n1,:)'); 
        corr_cy_RE_growth(i) =  corr(org_growth_RE_mat(i+2*n1,:)', org_growth_RE_mat(i+n1,:)'); 
        corr_Iy_RE_growth(i) =  corr(org_growth_RE_mat(i+4*n1,:)', org_growth_RE_mat(i+n1,:)');         
        corr_Pry_RE_growth(i) =  corr(org_growth_RE_mat(i+5*n1,:)', org_growth_RE_mat(i+n1,:)');   
        corr_PrH_RE_growth(i) =  corr(org_growth_RE_mat(i+5*n1,:)', org_growth_RE_mat(i,:)');         
     end
     
     ave_corr_Hy_growth = median(corr_Hy_growth); ave_corr_cy_growth = median(corr_cy_growth); 
     ave_corr_Iy_growth = median(corr_Iy_growth); ave_corr_Pry_growth = median(corr_Pry_growth); 
     ave_corr_PrH_growth = median(corr_PrH_growth); 
     
     ave_corr_Hy_RE_growth = median(corr_Hy_RE_growth); ave_corr_cy_RE_growth = median(corr_cy_RE_growth); 
     ave_corr_Iy_RE_growth = median(corr_Iy_RE_growth);ave_corr_Pry_RE_growth = median(corr_Pry_RE_growth)  ;
     ave_corr_PrH_RE_growth = median(corr_PrH_RE_growth); 
    display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')     
    display('correlations: growth rate data  ' )
    display('c,y I,y H,y    Pr,y Pr,H')
    display('Column 1: Learning model; Column 2 RE model')
    display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
     [ave_corr_cy_growth ave_corr_cy_RE_growth;
      ave_corr_Iy_growth ave_corr_Iy_RE_growth;         
      ave_corr_Hy_growth ave_corr_Hy_RE_growth;    
     ave_corr_Pry_growth ave_corr_Pry_RE_growth;
     ave_corr_PrH_growth  ave_corr_PrH_RE_growth] 
 
  std_corr_Hy_growth = std(corr_Hy_growth); std_corr_cy_growth = std(corr_cy_growth); 
     std_corr_Iy_growth = std(corr_Iy_growth); std_corr_Pry_growth = std(corr_Pry_growth); 
     std_corr_PrH_growth = std(corr_PrH_growth);
     
       display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')     
    display('standard deviation of correlations: growth rate data  ' )
    display('c,y I,y H,y    Pr,y Pr,H')
    display('Learning model')
    display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++') 
 [std_corr_cy_growth std_corr_Iy_growth std_corr_Hy_growth std_corr_Pry_growth std_corr_PrH_growth]
 
 % autocorr_org_growth = corr(org_growth_mat,0,2); 
 autocorr_c_growth = zeros(n1,1); 
 autocorr_I_growth = zeros(n1,1); 
 autocorr_Y_growth = zeros(n1,1); 
 autocorr_H_growth = zeros(n1,1); 
 autocorr_Pr_growth = zeros(n1,1); 
 
 autocorr_c_growth_RE = zeros(n1,1); 
 autocorr_I_growth_RE = zeros(n1,1); 
 autocorr_Y_growth_RE = zeros(n1,1); 
 autocorr_H_growth_RE = zeros(n1,1); 
 autocorr_Pr_growth_RE = zeros(n1,1); 
 
 autocorr_Y_growth_2lag =zeros(n1,1);
 autocorr_Y_growth_2lag_RE =zeros(n1,1);
 
for i=1:n1
    % orders H y c w inv Pr
    autocorr_c_growth(i)=corr(org_growth_mat(i+2*n1,1:end-1)', org_growth_mat(i+2*n1,2:end)');   
    autocorr_I_growth(i)=corr(org_growth_mat(i+4*n1,1:end-1)', org_growth_mat(i+4*n1,2:end)');   
    autocorr_Y_growth(i)=corr(org_growth_mat(i+n1,1:end-1)', org_growth_mat(i+n1,2:end)');   
    autocorr_H_growth(i)=corr(org_growth_mat(i,1:end-1)', org_growth_mat(i,2:end)');   
    autocorr_Pr_growth(i)=corr(org_growth_mat(i+5*n1,1:end-1)', org_growth_mat(i+5*n1,2:end)'); 
     
    autocorr_Y_growth_2lag(i)=corr(org_growth_mat(i+n1,1:end-2)', org_growth_mat(i+n1,3:end)');
    
    autocorr_c_growth_RE(i)=corr(org_growth_RE_mat(i+2*n1,1:end-1)', org_growth_RE_mat(i+2*n1,2:end)');   
    autocorr_I_growth_RE(i)=corr(org_growth_RE_mat(i+4*n1,1:end-1)', org_growth_RE_mat(i+4*n1,2:end)');   
    autocorr_Y_growth_RE(i)=corr(org_growth_RE_mat(i+n1,1:end-1)', org_growth_RE_mat(i+n1,2:end)');   
    autocorr_H_growth_RE(i)=corr(org_growth_RE_mat(i,1:end-1)', org_growth_RE_mat(i,2:end)');   
    autocorr_Pr_growth_RE(i)=corr(org_growth_RE_mat(i+5*n1,1:end-1)', org_growth_RE_mat(i+5*n1,2:end)');   
    
    autocorr_Y_growth_2lag_RE(i)=corr(org_growth_RE_mat(i+n1,1:end-2)', org_growth_RE_mat(i+n1,3:end)');
end

% Learning model
ave_autocorr_c_growth = median(autocorr_c_growth);
ave_autocorr_I_growth = median(autocorr_I_growth); 
ave_autocorr_Y_growth = median(autocorr_Y_growth); 
ave_autocorr_H_growth = median(autocorr_H_growth); 
ave_autocorr_Pr_growth = median(autocorr_Pr_growth); 

ave_autocorr_Y_growth_2lag = median(autocorr_Y_growth_2lag); 

% RE model
ave_autocorr_c_growth_RE = median(autocorr_c_growth_RE);
ave_autocorr_I_growth_RE = median(autocorr_I_growth_RE); 
ave_autocorr_Y_growth_RE = median(autocorr_Y_growth_RE); 
ave_autocorr_H_growth_RE = median(autocorr_H_growth_RE); 
ave_autocorr_Pr_growth_RE = median(autocorr_Pr_growth_RE); 

ave_autocorr_Y_growth_RE_2lag = median(autocorr_Y_growth_2lag_RE); 

display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
display('Autocorrelation (growth rates data)')
display('c, I, Y, H, Pr')
display('column 1 Learning model; column 2 RE model')
display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
[ave_autocorr_c_growth ave_autocorr_c_growth_RE;
 ave_autocorr_I_growth ave_autocorr_I_growth_RE; 
 ave_autocorr_Y_growth ave_autocorr_Y_growth_RE;
 ave_autocorr_H_growth ave_autocorr_H_growth_RE;
 ave_autocorr_Pr_growth ave_autocorr_Pr_growth_RE]
 
std_autocorr_c_growth = std(autocorr_c_growth);
std_autocorr_I_growth = std(autocorr_I_growth); 
std_autocorr_Y_growth = std(autocorr_Y_growth); 
std_autocorr_H_growth = std(autocorr_H_growth); 
std_autocorr_Pr_growth = std(autocorr_Pr_growth); 

display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
display('standard deviations of autocorrelation (growth rates data)')
display('c, I, Y, H, Pr')
display('Learning model')
display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')

[std_autocorr_c_growth std_autocorr_I_growth std_autocorr_Y_growth std_autocorr_H_growth std_autocorr_Pr_growth]

display('second order autocorrelation GDP growth learning vs RE')
[ave_autocorr_Y_growth_2lag ave_autocorr_Y_growth_RE_2lag]


sqrt(0.131*1.51^2)


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Kernel estimation of beliefs under learning
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% [fy,xi] = ksdensity(psiy(:,Tend)); 
% 
% %[fy_RE,xi] = ksdensity(median(org_growth_RE(n1+1:2*n1,:),2));
% 
% figure
% plot(xi,fy,'b-o'); 
% %hold on 
% %plot(xi, fy_RE, 'r-+')
% 
% 
% [fr,xi] = ksdensity(psir(:,Tend)); 
% figure
% plot(xi,fr); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the skewness
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



 skew_vec = skewness(var_level_detr');
 
  skew_H= median(skew_vec(1:n1));
  skew_y= median(skew_vec(n1+1:2*n1));
  skew_c= median(skew_vec(2*n1+1:3*n1));
  skew_w= median(skew_vec(3*n1+1:4*n1));
  skew_inv =median(skew_vec(4*n1+1:5*n1));
  skew_Pr =median(skew_vec(5*n1+1:6*n1));

  skew_vec_RE = skewness(var_level_detr_RE') ;
  
% orders H y c w inv Pr
  skew_H_RE= median(skew_vec_RE(1:n1));
  skew_y_RE= median(skew_vec_RE(n1+1:2*n1));
  skew_c_RE= median(skew_vec_RE(2*n1+1:3*n1));
  skew_w_RE= median(skew_vec_RE(3*n1+1:4*n1));
  skew_inv_RE =median(skew_vec_RE(4*n1+1:5*n1));
  skew_Pr_RE= median(skew_vec_RE(5*n1+1:6*n1));
  display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')    
  display('median skewness. Orders: H y c w inv Pr ')
  display('Column 1: Learning model; Column 2 RE model')
  display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')    
      
  [ skew_H  skew_H_RE;
     skew_y   skew_y_RE;
     skew_c  skew_c_RE;
     skew_w   skew_w_RE;
     skew_inv   skew_inv_RE;
     skew_Pr    skew_Pr_RE]
      

display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')    
display('HP filtered output volatility: Learning vs. RE')
display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')    

[median(STD_var_level_detr(n1+1:2*n1)) median(STD_var_level_detr_RE(n1+1:2*n1))]


xx= 0.0113*std_A/median(STD_var_level_detr(n1+1:2*n1));
yy = 0.0113*std_A/median(STD_var_level_detr_RE(n1+1:2*n1));


display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')  
display('Std of Productivity shock under Learning to match output volatility')
display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')  
 
xx
 
display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')  
display('Std of Productivity shock under RE to match output volatility')
display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')  

yy

display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
display('Percentage of productivity standard deviation reduced')
display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')

(xx-yy)/yy*100

display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
display('Output growth volatility under learning ')
display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')

y_growth_std400*xx/std_A


display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
display('Output growth volatility under RE ')
display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')

y_growth_std400_RE*yy/std_A

display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
display('')
display('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')

% 
% figure 
% 
% 
% subplot(2,1,1)
% plot(median(ForecastErrY(:,end-100:end),1))
% subplot(2,1,2)
% plot(median(ForecastErrRlevel1Q(:,end-100:end),1))

% figure 
% 
% subplot(2,1,1)
% set(gca,'fontsize',Font_size)
% plot(median(GDPgrowth(:,end-100:end),1),'r-o')
% hold on 
% plot(median(GDPgrowthfe(:,end-100:end),1))
% title('median GDP growth')
% legend('data','expectation')
% subplot(2,1,2)
% set(gca,'fontsize',Font_size)
% plot(median(interest(:,end-100:end),1),'r-o')
% hold on 
% plot(median(interestfe(:,end-100:end),1))
% title('median interest rate')
% legend('data','expectation')


display('Learning model performance: alternative gain')

display('gain ')

gainw

display('Amplification: output std & relative to RE')

(xx-yy)/yy*100

display('Amplification: I/Y std ratio & relative to RE')

   (ave_std_invyratio-ave_std_invyratio_RE)/ave_std_invyratio_RE
   
   
display('Amplification: H/Y std ratio & relative to RE')
(ave_std_Hyratio-ave_std_Hyratio_RE)/ave_std_Hyratio_RE 

display('Propagation Learning: GDP growth persistence')
ave_autocorr_Y_growth

% display('Forecast errors (rental rates) persistence')
% 
% ave_autocorr_FERlevel1Q

display('Relative volatility of the variance of productivity shocks to output growth volatility')

(1-0.34)^2* (xx*100)^2/(y_growth_std400*xx/std_A/4)^2


display('autocorrelation coefficient in 1-year ahead growth forecast errors')
median(autocorr_FEYlevel4Q)
display('standard deviation in autocorrelation coefficient in 1-year ahead growth forecast errors')
std(autocorr_FEYlevel4Q)

display('autocorrelation coefficient of forecast error in interest rate')
ave_autocorr_FERlevel1Q
display('standard deviation of autocorrelation coefficient in forecast error in interest rate')
std(autocorr_FERlevel1Q)

display('correlation between surprises to growth rates and revisions to potential growth forecasts')
median(corr_i1)

save coeff_mat_k_AL2000.mat coeff_mat_k_AL2000;
save coeff_mat_b_AL2000.mat coeff_mat_b_AL2000;
save coeff_mat_R_AL2000.mat coeff_mat_R_AL2000;
save coeff_mat_w_AL2000.mat coeff_mat_w_AL2000;

save y_tg2000.mat y_tg2000; 
save Dlogy2000.mat Dlogy2000; 
save Dlogw2000.mat Dlogw2000;
save DlogR2000.mat DlogR2000;
save Dlogk1999.mat Dlogk1999;
save Dlogb1999.mat Dlogb1999;
save Dlogk2000.mat Dlogk2000;
save Dlogb2000.mat Dlogb2000;
save  khat2000.mat  khat2000;
save  bhat2000.mat  bhat2000;

save what2000.mat what2000;
save RKhat2000.mat RKhat2000;

save yhat2000.mat yhat2000;
% save cab1996.mat cab1996;
% save cab1997.mat cab1997;
% save cab1998.mat cab1998;
% save cab1999.mat cab1999;
% save cab2000.mat cab2000;
cab_mat = cab; 

save cab_mat.mat cab_mat



