% ==================================================================
% This is a main program that calculates CS test results for   
% Korenok, Radchenko, Swanson "International Evidence on the Efficacy of
% new-Keynesian Models of Inflation Persistence"
% input is historical data and likelihood estimates of parameters
% output CS tables
% June 2006
% ==================================================================
%clc

clear;
% step 1. Load the data
% example US data; load raw data
Icomp = 1;
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';
end

rand_st =2051718; %sum(100*clock)*10 ;
randn('state',rand_st);
 
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;

outzaT = [];
outzaTstat = [];

    for I_subsample = 1:2; %cycle over samples
        for I_country = 1:13; % cycle over countries
            I_country
            outza_mes = [];
            outza_mesStat = [];
            for I_mes = 1:3; % cycle over measures
    
 [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;
    rpar = B(2:end,:);          % reduced form parameters
    
    % Solve SP, parameters: beta, lambda, v
    I_mod = 1;                  
    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);
    Cst_sp = Cst;
    [G1_sp, impact_sp, eu] = solveTeoreticalModel(spar, rpar, nvar, I_mod);
    
    % Solve SPI, parameters: beta, lambda, v
    I_mod = 2;                  
    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);    
    Cst_spi = Cst;
    [G1_spi, impact_spi, eu] = solveTeoreticalModel(spar, rpar, nvar, I_mod);

    % Solve SI, parameters: theta_3, xi, v
    I_mod = 3;
    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);
    Cst_si = Cst;
    [G1_si, impact_si, eu] = solveTeoreticalModel(spar, rpar, nvar, I_mod);

    indSC = [1;2;3]; %[1,2,3]-3 models; [1,2] - sp,spi; [1,3] - sp,si ; [3,2] - si,spi;
    aaSC = [50]; %[5;10;20;30;40;50]
    outza = CStest(Y(:,1),G1_sp,G1_spi, G1_si, impact_sp, impact_spi, impact_si, ...
                    Cst_sp, Cst_spi, Cst_si, indSC, aaSC);
    
    xlswrite(['SCres',num2str(I_mes),'.xls'], outza,I_country+(I_subsample-1)*13, 'C29:J29'); % save test results
    % summary table: 
    outza_mes = [outza_mes outza(end, end-2:end)];
    outza_mesStat = [outza_mesStat outza(end, 1:end-3)];
        end;
    outzaT = [outzaT; outza_mes];
    outzaTstat = [outzaTstat; outza_mesStat];
    end;
end;
xlswrite('CStablesF.xls', outzaT,3, 'C3:K28'); % save test results
xlswrite('CStablesF.xls', outzaTstat,4, 'C4:Q29'); % save SC statistics
load gong;
wavplay(y,Fs);