% Program that does calculations for "How Sticky Is Sticky Enough? 
% A Distributional and Impulse Response Analysis of New Keynesian DSGE
% Models"
% November 2004
% (c) Korenok and Swanson

clear;
% load raw data
% rgdp, cpi for 8 oecd countries, for list of countries see 'names'
load oecd_yp; 



% step 1. Define real data ld_pc_Y
ii = 21; % US
us_dat = oecd_dat(:,ii*2-1:ii*2);
us_dat = us_dat(us_dat(:,1)>0,:);
f_obs = 1960.25+(176-size(us_dat,1))*.25;
obs = (f_obs:.25:2004)';     % can-1961.25, che-1993.25, fra-1963.25, kor-1972.25, nze- 1961.5
rgdp = log(us_dat(2:end,1));
inf  = (us_dat(2:end,2) - us_dat(1:end-1,2))./us_dat(2:end,2);
rgdp_g = hpfilter(rgdp,1600)*rgdp; % output gap

% results: Data
ld_pc_Y(:,1) = rgdp_g(16:end) - mean(rgdp_g); %output gap
ld_pc_Y(:,2) = inf(16:end) - mean(inf); % inflation

% Summary statistic

% autocorrelation of inflation corr(pi_t, pi_t-1)
[z,c_inf] = corr(ld_pc_Y(:,2),ld_pc_Y(:,2),1); 
% autocorrelation of output corr(y_t, y_t-1)
[z,c_y] = corr(ld_pc_Y(:,1),ld_pc_Y(:,1),1); 
% cross correlation of inflation and output corr(y_t, pi_t)
[z,c_i_y] = corr(ld_pc_Y(:,1),ld_pc_Y(:,2),1);
cor_a = [c_inf(1); c_y(1); c_i_y(2)];

% cross correlation of inflation growth, pi(t+2)-pi(t-2), and output
[z,c_i_y4] = corr(ld_pc_Y(3:end-2,1),ld_pc_Y(5:end,2)-ld_pc_Y(1:end-4,2),1); 
% cross correlation of inflation growth, pi(t+4)-pi(t-4), and output
[z,c_i_y8] = corr(ld_pc_Y(5:end-4,1),ld_pc_Y(9:end,2)-ld_pc_Y(1:end-8,2),1); 
cros_cor_a = [c_i_y4(2);c_i_y8(2)];

% contingency table and statistic
[ct_a, ch2_a, pv_a ]= ct2(ld_pc_Y);
[ct_a_i, ch2_a_i, pv_a_i ]= ct2([ld_pc_Y(2:end,2) ld_pc_Y(1:end-1,2)]);
[ct_a_y, ch2_a_y, pv_a_y ]= ct2([ld_pc_Y(2:end,1) ld_pc_Y(1:end-1,1)]);

TT = size(ld_pc_Y,1); % define length of data sample


% step 2. Simulate the data according to alternative models

% Step 2.1. calibrated parameters:
a   = 2/3;  % labor share
b   = .99;  % beta, consumer discount factor
psi = 1;    % psi, intertemporal subst of labor 
s   = 1;    % intertemporal substitution of consumption
eps = 11;   % Dixit-Stiglitz parameter
th  = .5;  % price stickiness
p_a = 0;    % rho_a , ar parameter of technology process 
p_m = .5;   % rho_m , ar parameter exogenous money process
et  = 1;    % eta, parameter of money demand equation 
s_m = .007; % s.d. of residual for exogenous money process (MR)
s_a = .007; % s.d. of residual for technology process

% Step 2.2. intermediate calculations:
om = psi/a + 1/a -1;         % omega
psi_a = (om+1)/(s+om);

par = [a,b,psi,s,eps,th,p_a,p_m,et]';
Sig = [s_a 0; 0 s_m];

% Step 2.3. Models solution which later used for simulations
[G1_sp,impact_sp] = st_pr(par);  % sticky price model
[G1_spi,impact_spi] = st_pri(par);  % sticky price with indexation model
[G1_si,impact_si] = st_inf(par);  % sticky information model



%  Initial parameters: 

ls=[5;8;10;16;20];
n_ls = size(ls,1);
aa = [5;10;20;30;50;100];%[]
ind = [1,2,3]; % [1,2,3]-3 models; [1,2] - sp,spi; [1,3] - sp,si ; [3,2] - si,spi; 
n_aa = size(aa,1);
Unum=20;     % density of U grid 
bootnum=100;  % number of boot replications for CV calc 
simtot=1;

% for empirical application 

% grid size and coverage: 
Umin=min(ld_pc_Y)';
Umax=max(ld_pc_Y)';
Uinc=(Umax-Umin)./Unum;

outzaa=[];

% raw data assumed already to be in log differences
% define length of the simulated samples, given TT, where TT is actual data sample used 

outza = zeros(n_aa*n_ls,5+size(ind,2));
cor_s = zeros(n_aa,size(ind,2)*3);
cros_cor_s = zeros(n_aa,size(ind,2)*2);

%3. simulate data according to some different models, where SS is the length of the 
         
 % 3.1. simulated sample series
 x_sp = hist_sim(G1_sp,impact_sp,Sig,aa(end)*TT,[1,2,4]); %1-pi_t, 2-y_t, 4- a_t
 x_spi = hist_sim(G1_spi,impact_spi,Sig,aa(end)*TT,[1,2,4]);
 x_si = hist_sim(G1_si,impact_si,Sig,aa(end)*TT,[1,2,4]);

 % 3.2. simulated output gap and inflation
 YdatST =[x_sp(:,2)-psi_a*x_sp(:,3),x_spi(:,2)-psi_a*x_spi(:,3),x_si(:,2)-psi_a*x_si(:,3)];
 IdatST =[x_sp(:,1),x_spi(:,1),x_si(:,1)]; % inflation
 YdatST = YdatST(:,ind);
 IdatST = IdatST(:,ind);
 Oboy = zeros(n_aa*n_ls,size(YdatST,2)-1);

for tsi = 1:n_aa;
   SS = aa(tsi)*TT;
   YdatS = YdatST(1:SS,:);
   IdatS = IdatST(1:SS,:);
         % 3.3. additional statistics:
         
         % autocorrelation of inflation
         [z,c_inf_s] = corr(IdatS,IdatS,2); 
         % autocorrelation of output
         [z,c_y_s] = corr(YdatS,YdatS,2); 
         % cross correlation of inflation and output
         [z,c_i_y_s] = corr(IdatS,YdatS,2); 
         
         cor_s(tsi,:) = [c_inf_s(1,:) c_y_s(1,:) c_i_y_s(2,:)];
         
         [z,c_i_y4_s] = corr(YdatS(3:end-2,:),IdatS(5:end,:)-IdatS(1:end-4,:),2); 
         % cross correlation of inflation growth, pi(t+2)-pi(t-2), and output
         [z,c_i_y8_s] = corr(YdatS(5:end-4,:),IdatS(9:end,:)-IdatS(1:end-8,:),2); 
         % cross correlation of inflation growth, pi(t+4)-pi(t-4), and output
         cros_cor_s(tsi,:) = [c_i_y4_s(2,:) c_i_y8_s(2,:)];

         % 3.4. now, caculate all stats over U range 
         [outttu simF] = zee2o(ld_pc_Y(:,1),ld_pc_Y(:,2),YdatS,IdatS,Umin,Uinc,Unum);
         % average up across the u, rational for averagin described on p.6 
         oboy=mean(outttu)';
         % 3.5. calculate statistics 
         ZZ=ZTstats(oboy);
         Oboy((tsi-1)*n_ls+1:tsi*n_ls,:) = repmat(oboy',5,1); 
         outza((tsi-1)*n_ls+1:tsi*n_ls,[1,6:6+size(IdatST,2)-1])=[ZZ*ones(n_ls,1) repmat(mean(simF),n_ls,1)]; 
end;

%        4.  Now, bootstrap the critical values
ZZZb1=zeros(bootnum,n_aa*n_ls); ZZZb2=zeros(bootnum,n_aa*n_ls);

for booti = 1:bootnum;
    % 4.1. simulated sample series
    x_sp = hist_sim(G1_sp,impact_sp,Sig,aa(end)*TT,[1,2,4]); %1 - pi_t, 2-y_t, 4- a_t
    x_spi = hist_sim(G1_spi,impact_spi,Sig,aa(end)*TT,[1,2,4]);
    x_si = hist_sim(G1_si,impact_si,Sig,aa(end)*TT,[1,2,4]);
         
    % 4.2. simulated output gap
    YdatS_bT =[x_sp(:,2)-psi_a*x_sp(:,3),x_spi(:,2)-psi_a*x_spi(:,3),x_si(:,2)-psi_a*x_si(:,3)];
    IdatS_bT =[x_sp(:,1),x_spi(:,1),x_si(:,1)];
    YdatS_bT = YdatS_bT(:,ind);
    IdatS_bT = IdatS_bT(:,ind); 
    u = rand(aa(end)*TT,1);  

    for tsi = 1:n_aa;
        SS = aa(tsi)*TT;
        YdatS_b = YdatS_bT(1:SS,:);
        IdatS_b = IdatS_bT(1:SS,:);

        for li = 1:n_ls;
            l_val=ls(li); % l_val for actual data 

            % now, make l_vals for simulated data
            l_valS = aa(tsi)*ls(li);
            % assume that the "l_val" values are for the actual data -- for the simulated, mutiply them
            % by the factor 1, 2, or 5, as SS is made by the same factors multiplying TT 

            bdatAct = bootk1((ld_pc_Y),l_val,u);
            bdatSim = bootk1([YdatS_b IdatS_b],l_valS,u);

            % 4.3.  I carry out the usual boot approach 

            outttu = zee2o(bdatAct(:,1),bdatAct(:,2),bdatSim(:,1:size(YdatS_b,2)),...
                     bdatSim(:,size(YdatS_b,2)+1:end),Umin,Uinc,Unum);
            boboy1=mean(outttu)-Oboy(li+(tsi-1)*n_ls,:); 
            % these are the usual boot stats with the booted minus the actual 

            % 4.4. II carry out the alternative boot for S growing faster than T 

            outttu = zee2o(bdatAct(:,1),bdatAct(:,2),YdatS_b,IdatS_b,Umin,Uinc,Unum);
            boboy2=mean(outttu)-Oboy(li+(tsi-1)*n_ls,:); 
         
            % these are the alt boot stats with the booted minus the actual, ut where
            % only the actual data are resampled, and not the booted, i.e. S grow faster T 
             % calculate statistics 

           ZZb1=ZTstats(boboy1');
           ZZZb1(booti,li+(tsi-1)*n_ls)=ZZb1;

           ZZb2=ZTstats(boboy2');
           ZZZb2(booti,li+(tsi-1)*n_ls)=ZZb2;            
        end
    end
    if rem(booti,10)==0; booti 
    end;
end  
         
% 4.5. construct critical values 

temp=sort(ZZZb1,1);
cv1_10=temp(floor(0.90*size(ZZZb1,1)),:);
cv1_5=temp(floor(0.95*size(ZZZb1,1)),:);

temp=sort(ZZZb2,1);
cv2_10=temp(floor(0.90*size(ZZZb2,1)),:);
cv2_5 =temp(floor(0.95*size(ZZZb2,1)),:);


         % 5. Results: 
         % 5.1. statistics and critical values
         
         outza(:,2:5)=[cv1_5' cv1_10' cv2_5' cv2_10'];
         
         % Table: Distributional Accuracy Tests Based on the Comparison of Historical and Simulated  
         % inflation and output gap
         display('Distributional Accuracy Tests')
         outza  
         
         % 5.2. Quantiles, actual vs simulated:
         p = [.1 .2 .3 .4 .5 .6 .7 .8 .9];
         
         q_a = quantile(ld_pc_Y, p); % actual data
         q_ys = quantile(YdatS, p); % simulated output
         q_is = quantile(IdatS, p); % simulated inflation
         
         % Historical and Simulated Distribution Quantiles for Inflation and the Output gap
         display('Historical and Simulated Distribution Quantiles')
         display('output gap')
         q_Y  = [q_a(:,1)'; q_ys'] % output table
         display('inflation')
         q_I  = [q_a(:,2)'; q_is'] % inflation table
         
         % 5.3. Contingency tables
         % out_g, inf
         
         [ct_sp, ch2_sp, pv_sp ]= ct2([YdatS(:,1),IdatS(:,1)]);
         [ct_spi, ch2_spi, pv_spi]= ct2([YdatS(:,2),IdatS(:,2)]);
         [ct_si, ch2_si, pv_si ]= ct2([YdatS(:,3),IdatS(:,3)]);
         
         
         % inf_t, inf_t-1
         [ct_sp, ch2_sp, pv_sp ]= ct2([IdatS(2:end,1),IdatS(1:end-1,1)]);
         [ct_spi, ch2_spi, pv_spi]= ct2([IdatS(2:end,1),IdatS(1:end-1,2)]);
         [ct_si, ch2_si, pv_si ]= ct2([IdatS(2:end,1),IdatS(1:end-1,3)]);
         
         % y_t, y_t-1
         [ct_sp, ch2_sp, pv_sp ]= ct2([YdatS(2:end,1),YdatS(1:end-1,1)]);
         [ct_spi, ch2_spi, pv_spi]= ct2([YdatS(2:end,1),YdatS(1:end-1,2)]);
         [ct_si, ch2_si, pv_si ]= ct2([YdatS(2:end,1),YdatS(1:end-1,3)]);
         
         % Table: Directional Accuracy of Simulated Ination and the Output Gap Paths
         % Inflaiton vs Output gap
         display('Directional Accuracy')
         display('Inflation, Output')
         CT_yi = [[ct_a(:);pv_a],[ct_sp(:);pv_sp],[ct_spi(:);pv_spi]...  
               [ct_si(:);pv_si]]
         display('Inflation vs Lagged inflation')
         CT_ii = [[ct_a_i(:);pv_a_i],[ct_sp(:);pv_sp],[ct_spi(:);pv_spi]...  
               [ct_si(:);pv_si]]
         display('Output gap vs Lagged output gap')
         CT_yy = [[ct_a_y(:);pv_a_y],[ct_sp(:);pv_sp],[ct_spi(:);pv_spi]...  
               [ct_si(:);pv_si]]
           
         %5.4
         % Table: Autocorrelation and Cross Correlation: Inflation and the Output Gap
         display('Autocorrelation and Cross Correlations')
         cor_s
         
         % 5.5
         % Table: Correlation Between Output Gap and Change in Inflation: Acceleration Phenomena
         display('Correlation Between Output Gap and Change in Inflation: Acceleration Phenomena')
         cros_cor_s
         
         % 5.6. Actual vs simulated densities:
         [x_a,den_a]=kden(ld_pc_Y); % actual density
         obs_kd = (1:10:size(YdatS,1))';  
         [x_ys,den_ys]=kden(YdatS(obs_kd,:)); % simulated output density
         [x_is,den_is]=kden(IdatS(obs_kd,:)); % simulated inflation density
         % plot
         figure; plot(x_a(:,1),den_a(:,1),x_ys(:,1),den_ys(:,1),x_ys(:,2),den_ys(:,2),...
             x_ys(:,3),den_ys(:,3)); legend('act','sp','spi','si'); 
          
         figure; plot(x_a(:,2),den_a(:,2),x_is(:,1),den_is(:,1),x_is(:,2),den_is(:,2),...
             x_is(:,3),den_is(:,3)); legend('act','sp','spi','si');