% ==================================================================
% This program computes residuals autocorrelation, RMSE, MAE, regress 
% the constant model residuals on SP, SPI and SI
% Korenok, Radchenko, Swanson "International Evidence on the Efficacy of
% new-Keynesian Models of Inflation Persistence"
% 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\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
Table= zeros(nC*4,12,3); 

    Res_count1 = [];
    Res_count2 = [];
    Res_count3 = [];        
    Res_count4 = [];     
    for I_subsample = 1:2; %cycle over samples   
        for I_country = 1:13; % cycle over countries
            Mod1 = [];
            Mod2 = [];
            Mod3 = [];
            Mod4 = [];
            for I_mes = 1:3; % cycle over measures
                vtm = [];
                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);
    % calculate fitted values
    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);
    vtm = [vtm Y(:,1)-vt(1,:)'];
                end;
                ind = (2:10:size(Y,1)-1)';
                Yhat_var = X*B'; %VAR forecast
                dates1 = dates(end-size(Y,1)+2:end);
                vtm = ((1+mX(:,1)+ vtm(2:end,:)).^4 - 1)*100;
                %Yhat_var_pi = ((1+mX(:,1)+Yhat_var(2:end,1)).^4-1)*100;% VAR forecast
                Yhat_var_pi = ((1+mX(:,1)*ones(size(Yhat_var(2:end,1),1),1)).^4-1)*100;% constant inf 
                Y_pi = ((1+mX(:,1) + Y(2:end,1)).^4 - 1)*100;
%                figure; plot(dates1,Y_pi,'-k',dates1,Yhat_var_pi,'+',...
%                     dates1, vtm(:,1),'x',dates1, vtm(:,2),':', dates1, vtm(:,3),'--'...
%                     ,'LineWidth',1); 
%                legend('Historical','VAR','Sticky Price','Sticky Price with Indexation',...
%                    'Sticky Information');  axis([dates1(1) dates1(end) 0 13 ])
                
                ehat = repmat(Y_pi,1,4) - [Yhat_var_pi, vtm];
                Yhat = [Yhat_var_pi, vtm];
                acf_ac   = [];              % sample autocorrelations
                acf_Q    = [];              % sample q-stat
                acf_p    = [];              % p-values for Q-stat
                cor_epi  = [];
                res_mes = [];
                for ik = 1:4;
                    SACF = acf(ehat(:,ik), 1);      
                    acf_ac = [acf_ac SACF.ac'];
                    acf_Q  = [acf_Q SACF.qstat'];
                    acf_p  = [acf_p SACF.prob'];
                    e_mean = Y_pi-mX(:,1)*ones(size(Y_pi,1),1);
                    cor_temp = corrcoef([e_mean ehat(:,ik)]);
                    %cor_temp = corrcoef([Y_pi Yhat(:,ik)]);
                    cor_temp = cor_temp(2,1);
                    cor_epi  = [cor_epi cor_temp];
                end
                ehat_m = [mean(ehat)]';                 % mean error
                RMSE   = [sqrt(mean(ehat.^2))];        
                MAE    = [mean(abs(ehat))];
                RegrRes  = [];
                for i = 1:3;
                    X1 = [ones(size(ehat,1),1) ehat(:,i+1)];
                    Y1 = Y_pi;
                    B1       = X1\Y1;                      
                    err1     = Y1 - X1*B1;
                    R2 = 1 - (err1'*err1)/(Y1'*Y1);
                    RegrRes = [RegrRes R2];
                end
    
                    
                Mod1   = [Mod1 RMSE];   
                Mod2   = [Mod2 MAE];    
                Mod3   = [Mod3 cor_epi];
                Mod4   = [Mod4 RegrRes];
%                Mod1   = [Mod1 ehat_m(1) RMSE(1) acf_ac(1)];   %VAR
%                Mod2   = [Mod2 ehat_m(2) RMSE(2) acf_ac(2)];   %SP
%                Mod3   = [Mod3 ehat_m(3) RMSE(3) acf_ac(3)];   %SPI
%                Mod4   = [Mod4 ehat_m(4) RMSE(4) acf_ac(4)];   %SI                

                igraph = 0;
                if igraph ==1;
                figure; plot(dates1,ehat(:,1),'+',...
                     dates1, ehat(:,2),'x',dates1,ehat(:,3),':', dates1, ehat(:,4),'--'...
                     ,'LineWidth',1); 
                legend('VAR','Sticky Price','Sticky Price with Indexation',...
                    'Sticky Information'); axis([dates1(1) dates1(end) -5 9 ])
                end
    
            end;
        Res_count1 = [Res_count1; Mod1(:,1:4) Mod2(:,1:4) Mod3(:,1:4)];     % output gap
        Res_count2 = [Res_count2; Mod1(:,5:8) Mod2(:,5:8) Mod3(:,5:8)];     % labor share
        Res_count3 = [Res_count3; Mod1(:,9:12) Mod2(:,9:12) Mod3(:,9:12)];  % unemployment
        Res_count4 = [Res_count4; Mod4];
        %        Res_count4 = [Res_count4; Mod4];
    end;
end;
ResCorr = [Res_count1(:,end-2:end) Res_count2(:,end-2:end) Res_count3(:,end-2:end)];
ResOLS = [Res_count1(:,end-2:end) Res_count2(:,end-2:end) Res_count3(:,end-2:end)];

%xlswrite('res_ehat2', Res_count1, 1, 'C3:N30');
%xlswrite('res_ehat2', Res_count2, 2, 'C3:N30');
%xlswrite('res_ehat2', Res_count3, 3, 'C3:N30');
%xlswrite('res_ehat2', ResCorr , 1, 'C3:K28');
%xlswrite('res_ehat2', ResCorr , 2, 'C3:K28');
%xlswrite('res_ehat2', Res_count4 , 4, 'C3:K28');