% ==================================================================
% This program produces MLE parameter estimates (beta constrained at .99, 
% and 1 lag is used for all models/countries/measures) 
% 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 = 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
%indCount = [1,4,10,11,12,15,18,22,23,24,28,30]; % without Italy
%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);



% ======================================
%   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); 

for I_mod = 2:2; % cycle over models
    res_sub = [];
    for I_subsample = 1:2; %cycle over samples
        for I_country = 13:13; % cycle over countries
            res_mes = [];
            for I_mes = 2:2; % cycle over measures
                
    
 X  = 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)    
    spar = [b; la; th; xi];     % structural parameters
    rpar = B(2:end,:);          % reduced form parameters

    [G1, impact, eu] = solveTeoreticalModel(spar, rpar, nvar, I_mod);
    
    if I_mod == 1;                  % Parameters for SP: beta, lambda, v
        spar0   = [.001; .003]; % Starting values
        if I_country == 11; spar0 = [.005;.003];end;
        [G1, impact, eu] = solveTeoreticalModel([spar(1); spar0(1); spar(3:4)], rpar, nvar, I_mod);
        if sum(eu)~=2; spar0(1)=-spar0(1);end

        lb = [-.3; .0001];        % non-binding lower bound 
        ub = [ .5; 1];      % non-binding upper bound    
    elseif I_mod == 2;                  % Parameters for SPI: beta, lambda, v
        spar0   = [.0001; .003];   % Starting values
        [G1, impact, eu] = solveTeoreticalModel([spar(1); spar0(1); spar(3:4)], rpar, nvar, I_mod);
        if sum(eu)~=2; spar0(1)=-spar0(1);end

        lb = [ -.3; .0001];       %non-binding lower bound 
        ub = [.5; 1];        %non-binding upper bound    
    else
        %   Parameters for SI: theta, xi, v
        spar0   = [.5; .001; .009];  % Starting values
        [G1, impact, eu] = solveTeoreticalModel([spar(1:2) spar0(1:2)], rpar, nvar, I_mod);
        if sum(eu)~=2; spar0(2)=-spar0(2);end
        lb = [.01; -.5; .001];       %non-binding lower bound 
        ub = [1; .5; .5]; %non-binding upper bound 
        if I_country > 12; lb(2) = -1; ub(2) = 1; end;  
    
    end
        [Par_est, Hes, lik] = LikModel2(spar0, Y, X, Cst, rpar, B, nvar,...
                        I_mod, lb, ub, spar);    
        Par_se      = sqrt(diag(inv(Hes)));
        res_par     = [Par_est' -lik];
        res_se      = [Par_se' 0 ];
        res_all     = [res_par; res_se];
        res_mes     = [res_mes res_all];
              end;
              res_sub = [res_sub; res_mes];
        end;
    end;
    if I_mod<3;
    xlswrite('res_mle2', res_sub,I_mod, 'C3:K54');
    xlswrite('res_mle2', lb,I_mod, 'N3:N4');
    xlswrite('res_mle2', ub,I_mod, 'O3:O4');
    xlswrite('res_mle2', spar0,I_mod, 'N7:N8');
    else;
    xlswrite('res_mle2', res_sub,I_mod, 'C3:N54');
    xlswrite('res_mle2', lb,I_mod, 'Q3:Q5');
    xlswrite('res_mle2', ub,I_mod, 'R3:R5');
    xlswrite('res_mle2', spar0,I_mod, 'Q8:Q10');
    end
    %res_sub
end;