clear; clc; close all;
cd('/Users/julioortiz/Dropbox/ME/replication/');
addpath '/Users/julioortiz/Dropbox/ME/replication/code/functions';

macrovar=["rgdp" "pgdp" "cpi" "emp" "unemp" "rconsum" "rnresin" "rresinv" "rfedgov" "rslgov" "housing" "indprod" "tbill" "tbond"];
LRtest = zeros(length(macrovar),3,2);

for i=1:length(macrovar)
    
clear fe fr lfe Y
%% First, solve actuals
dat=readmatrix(strcat('processed_data/',macrovar(i),'_actuals.csv'));
year=dat(:,1);
quarter=dat(:,2);
act_rt=dat(:,3);
tStart=tic;
[xvec,fval]=filterfnc_ar2_H(act_rt);
tEnd=toc(tStart);
fprintf('Parameters estimated in %d minutes and %f seconds.\n',floor(tEnd/60), rem(tEnd,60));

%% Next, solve panel MLE from 1992-2019
dat=readmatrix(strcat('processed_data/spf9219_',macrovar(i),'_spell8_fefrlfe.csv'));
id = dat(:,1);
year = dat(:,2);
quarter = dat(:,3);
fe = dat(:,4);
fr = dat(:,5);
lfe = dat(:,7);
fcast0q=dat(:,8);
fcast1q=dat(:,9);

if macrovar(i)=="cpi"
    first = 1994;
elseif macrovar(i)=="emp"
    first = 2004;
else
    first = 1992;
end
last = 2019;

if macrovar(i)=="tbond"
    T = (last-first+1)*4-1;
elseif macrovar(i)=="cpi"
    T = (last-first)*4+2;
else
    T = (last-first+1)*4;
end

N = size(dat(:,1),1)/T;

fe=reshape(fe,[T,N])';
fr=reshape(fr,[T,N])';
lfe=reshape(lfe,[T,N])';

Y(:,:,1)=fe;
Y(:,:,2)=fr;
Y(:,:,3)=lfe;

tStart=tic;
[xvec_RE,fval_RE]=filtervecfnc_RE(Y,N,T,xvec(2:4)); %RE model
[xvec_OC,fval_OC]=filtervecfnc_OC(Y,N,T,xvec(2:4),xvec_RE); %OC model
[xvec_DE,fval_DE]=filtervecfnc_DE(Y,N,T,xvec(2:4),xvec_RE); %DE model
[xvec_ME,fval_ME]=filtervecfnc_ME(Y,N,T,xvec(2:4),xvec_RE); %ME model
tEnd=toc(tStart);
fprintf('Parameters estimated in %d minutes and %f seconds.\n',floor(tEnd/60), rem(tEnd,60));

%% Now do Vuong (1989) non-nested test

% ME vs. RE
[omegatilde,LL_RE,LL_ME,term1,LL_RE_cum,LL_ME_cum]=LL_vuong_MERE(Y,N,T,xvec,xvec_RE,xvec_ME);
Vstat = (LL_ME_cum - LL_RE_cum)*term1^(-0.5)/(omegatilde^0.5);
LRtest(i,1,1) = Vstat;
LRtest(i,1,2) = (1-normcdf(Vstat));

% ME vs. OC
[omegatilde,LL_OC,LL_ME,term1,LL_OC_cum,LL_ME_cum]=LL_vuong_MEOC(Y,N,T,xvec,xvec_RE,xvec_OC,xvec_ME);
Vstat = (LL_ME_cum - LL_OC_cum)*term1^(-0.5)/(omegatilde^0.5);
LRtest(i,2,1) = Vstat;
LRtest(i,2,2) = (1-normcdf(Vstat));

% ME vs. DE
[omegatilde,LL_DE,LL_ME,term1,LL_DE_cum,LL_ME_cum]=LL_vuong_MEDE(Y,N,T,xvec,xvec_RE,xvec_DE,xvec_ME);
Vstat = (LL_ME_cum - LL_DE_cum)*term1^(-0.5)/(omegatilde^0.5);
LRtest(i,3,1) = Vstat;
LRtest(i,3,2) = (1-normcdf(Vstat));

end 
LRstat = LRtest(:,:,1);
LRpval = LRtest(:,:,2);
[macrovar' LRstat]