% ==================================================================
% This program computes theoretical ACF 
% Korenok, Radchenko, Swanson "International Evidence on the Efficacy of
% new-Keynesian Models of Inflation Persistence"
% June, 2006
% ==================================================================

clear;
% step 1. Load the data
% example US data; load raw data
Icomp = 0;
if Icomp == 1;
    dataplace = 'C:\Documents and Settings\Faculty\Desktop\current projects\infoIV\data';
    progplace = 'C:\Documents and Settings\Faculty\Desktop\current projects\infoIV\mlab';
else
    dataplace = 'H:\My Documents\SPSImodels\data';
    progplace = 'H:\My Documents\SPSImodels\matlab\MLE';
end

 
indCount = [1,4,10,11,12,15,17,18,22,23,24,28,30]; % selecting countries
%Aust, Can, Fin, Fra, UK, Ire, Ita, Jap, Neth, Nor, NZ, Swe, US
cd(dataplace);
filename = 'OECDEconomicOutlook77.xls';     range   = 'B7:AE188'; 
rgdp    = xlsread(filename, 2, range); rgdp = rgdp(:,indCount);
yg      = xlsread(filename, 5, range); yg = yg(:,indCount);
pop     = xlsread(filename, 8, range); pop = pop(:,indCount);
gdpDef  = xlsread(filename, 9, range); gdpDef = gdpDef(:,indCount);
ulc     = xlsread(filename, 3, range); ulc = ulc(:,indCount);
unemp   = xlsread(filename, 6, range); unemp = unemp(:,indCount);
intr    = xlsread(filename, 12, range);intr = intr(:,indCount);
cd(progplace);

par_sp = xlsread('res_mle2', 1, 'C3:K54');
par_spi = xlsread('res_mle2', 2, 'C3:K54');
par_si  = xlsread('res_mle2', 3, 'C3:N54');

% ======================================
%   Variables transformations
% ======================================
dates = (1960.25:.25:2005.5)';
[tC, nC] = size(yg);
rulc= zeros(tC,nC); unempA = rulc;
intrQuart = rulc; intrA = rulc; 
piDefQuart =zeros(tC-1,nC); piDefYear = zeros(tC-4,nC);
for i=1:nC
    yg(yg(:,i)==-10000,i)    = 0;                   % this is just for output gap, excel related
    yg(yg(:,i)<-10,i)        = yg(yg(:,i)<-10,i)+100000;
    yg(:,i)                  = yg(:,i)/100;
    ind  = find(ulc(:,i)~=-10000);
    rulc(ind,i) = log(ulc(ind,i)) - log(gdpDef(ind,i)); % real unit labor cost (labor share) 
    %dulc(ind,i) = log(ulc(ind(2:end),i)) - log(ulc(ind(1:end-1),i)); % quarterly inflation in ulc
    ind  = find(unemp(:,i)~=-10000);
    unempA(ind,i) = unemp(ind,i)/100;
    unemp(:,i) = unempA(:,i);
    ind  = find(intr(:,i)~=-10000);
    intrQuart(ind,i)= (1+intr(ind,i)/100).^(1/4)-1;     % quarterly interest rate
    intrA(ind,i) = intr(ind,i)/100;
    intr(:,i) = intrA(:,i);
    %figure; plot(dates,intr(:,i));
    ind  = find(gdpDef(:,i)~=-10000);
    [piDefQuart(ind(1:end-1),i) piDefYear(ind(1:end-4),i)] = Quart_trans(gdpDef(ind,i));       % quaterly inflation
    ind  = find(rgdp(:,i)~=-10000);
    [rgdpPerCapitaGrQuart(ind(1:end-1),i) rgdpPerCapitaGrYear(ind(1:end-4),i)]...
        = Quart_trans((rgdp(ind,i)./pop(ind,i))); % quarterly per capita output
    %figure; plot(dates(5:end),rgdpPerCapitaGrYear(:,i));
end

% ===========================================
%   Controls for estimation
% ===========================================
I_mes       =1;        % output gap measure: 
                        % 1  - "output" gap; 2 - unit labor cost; 3- unemployment
I_mod       = 1;        % 1- sp; 2 - spi;   3 - si                        
I_freq      = 1;        % 1- quarterly defenition; 2 - annual defenition
I_subsample = 2;        % 1 - full sample
I_country   = 13;

% Table structure: top to bottom: countries full sample, countries second
% sample; left to rigth parameter estimates for ouput gap measure,
% parameter estimates for labor share measure, parameter estimates for
% negative unemployment measure; one table for sp, one table for spi, one
% table for si
acfTable =[];                   acf_e_Table =[];
t1Table  =[];                   t2Table =[];
t3Table  =[];


    for I_subsample = 1:2; %cycle over samples
        for I_country = 13:13; % cycle over countries
            acfM = [];            rho_eM = [];
            t1M  = [];            t2M  = [];                t3M = [];

            for I_mes = 1:3; % cycle over measures
                acf_y   = [];           acf_e   = [];               
                t1_a    = [];           t2_a    = [];       t3_a = [];
                for I_mod = 1:3; % cycle over models
    
 [X mX] = def_X(piDefQuart,  yg, rulc, unemp, rgdpPerCapitaGrQuart, intrQuart,...
 piDefYear, rgdpPerCapitaGrYear, intr,dates, I_freq, I_subsample, I_mes, I_country);
    % ===========================================
    %       VAR estimation
    % ===========================================
    nlag    = 1;                % number of lags in VAR
    [Y, X]  = var_data(X, nlag);
    [T,nvar]= size(Y);
    B       = X\Y;                      
    err     = Y - X*B;
    Sig_var = err'*err/T;
    C       = chol(cov(err))';    % Y[t] = X[t]*B + w[t]*C
                        % Calculating Cholesky decompos. for struct. model
    Cs      = chol(cov(err(:,2:end)))';     % covariance of tw[t]
    se_v    = C(1,1);                       % covariance of v[t]
    Cst     = eye(nvar); 
    Cst(1,1)= se_v;             Cst(2:end,2:end) = Cs;
    B = [B'; eye(nvar*(nlag-1)) zeros(nvar*(nlag-1),nvar)]; % matrix in companion form

    % ===============================================
    %  parameters in spar are used as reference values
    % ===============================================
    [b, la, th, xi] = def_spar; % Starting values (prior expectation)    
    ind_r = 2*(I_country+(I_subsample-1)*13)-1;
    if I_mod == 1;                  % Parameters for SP: beta, lambda, v
        spar = [b; par_sp(ind_r,I_mes*3-2); th; xi];     % structural parameters
        Cst(1,1) = par_sp(ind_r,I_mes*3-1);
    elseif I_mod == 2;                  % Parameters for SPI: beta, lambda, v
        spar = [b; par_spi(ind_r,I_mes*3-2); th; xi];     % structural parameters
        Cst(1,1) = par_spi(ind_r,I_mes*3-1);    
    else
        spar = [b; la; par_si(ind_r,I_mes*4-3); par_si(ind_r,I_mes*4-2)];     % structural parameters
        Cst(1,1) = par_si(ind_r,I_mes*4-1);
    end

    
    
    rpar = B(2:end,:);          % reduced form parameters
    % solve the model
    [G1, impact, eu] = solveTeoreticalModel(spar, rpar, nvar, I_mod);
    
    H    = [eye(nvar) zeros(nvar,size(G1,1)-nvar)]; %"selection" matrix
    xtt1 = zeros(size(G1,1),1); 
    K_x = size(X,2);            
    xtt1(1:K_x) = (X(1,:)*B')';
    [f vt]     = LogLik_SP1(Y, H , G1, impact, Cst*Cst',xtt1);
    
    rho_e  = acf(vt(1,:)', 1);     
    rho_e = rho_e.ac;        % residuals acf
    
    Rho = compACF(G1, impact, Cst, Y, 1);
    t1  = Rho(1,1)*sqrt(T);             % should be N(0,1)
    t2  = T*(T+2)*(Rho(1,1)^2)/(T-1);   % Q-stat for only one autocorrelation
    t3  = T*(T+2)*(rho_e(1,1)^2)/(T-1);   % Q-stat for only one autocorrelation
    
    acf_y   = [acf_y Rho(1,1)];         acf_e   = [acf_e rho_e];
    t1_a    = [t1_a t1];                t2_a    = [t2_a t2];
    t3_a    = [t3_a t3];
          end;
    acfM = [acfM acf_y];            rho_eM  = [rho_eM acf_e];
    t1M  = [t1M t1_a];              t2M     = [t2M t2_a];
    t3M  = [t3M t3_a];
          end;
          acf_var = acf(err(:,1),1); acf_var = acf_var.ac;
          
    acfM = [Rho(2,1) acfM];         rho_eM  = [acf_var rho_eM];
    t1M  = [Rho(2,1)*sqrt(T) t1M];  t2M     = [T*(T+2)*(Rho(2,1)^2)/(T-1) t2M];
                                    t3M     = [T*(T+2)*(acf_var(1,1)^2)/(T-1) t3M];
                                    
         acfTable = [acfTable; acfM];       acf_e_Table = [acf_e_Table; rho_eM];
         t1Table  = [t1Table; t1M];         t2Table     = [t2Table; t2M];         
         t3Table  = [t3Table; t3M];  
    end;
end;
xlswrite('res_ehat2', acfTable , 7, 'C3:L28');
xlswrite('res_ehat2', t1Table , 8, 'C3:L28');
xlswrite('res_ehat2', t2Table , 9, 'C3:L28');
xlswrite('res_ehat2', acf_e_Table , 10, 'C3:L28');
xlswrite('res_ehat2', t3Table , 11, 'C3:L28');