clear; clc; close all;

cd('/Users/julioortiz/Dropbox/ME/replication/');

macrovar=["rgdp" "pgdp" "cpi" "emp" "unemp" "rconsum" "rfedgov" "rslgov" "rnresin" "rresinv" "tbill" "tbond" "housing" "indprod"];

Encompassing_Weights=zeros(length(macrovar),4);
fund_parameter_estimates = zeros(length(macrovar),4);
rhohatvec = zeros(length(macrovar),1);

for i=1:length(macrovar)

%% First, solve for fundamental parameters from 1992-2005
clear year quarter act_rt Y Y2 dat dat2 fe fr femn N T N2 T2 id year quarter

dat=readmatrix(strcat('processed_data/',macrovar(i),'_actuals.csv'));
year=dat(:,1);
quarter=dat(:,2);
act_rt=dat(:,3);

if macrovar(i)=="cpi"
    Y=act_rt(1:51); %1994-2006
elseif macrovar(i) == "emp"
    Y=act_rt(49:80); %2004-2011
else
    Y=act_rt(1:56); %1992-2005
end

tStart=tic;
[xvec,fval]=filterfnc_ar2_H(Y);
fund_parameter_estimates(i,:) = xvec;
tEnd=toc(tStart);
fprintf('Parameters estimated in %d minutes and %f seconds.\n',floor(tEnd/60), rem(tEnd,60));

%% First, solve via MLE from 1992-2005
dat=readmatrix(strcat('processed_data/spf9219_',macrovar(i),'_oos1.csv'));
id = dat(:,1);
year = dat(:,2);
quarter = dat(:,3);
fe = dat(:,4);
fr = dat(:,5);
lfe = dat(:,6);
f0data = dat(:,7);
f1data = dat(:,8);

if macrovar(i)=="cpi"
    first = 1995;
    last = 2006;
    T = (last-first+1)*4+2;
elseif macrovar(i) == "emp"
    first=2004;
    last=2011;
    T = (last-first+1)*4;
elseif macrovar(i)=="tbond"
    first = 1992;
    last = 2005;
    T = (last-first+1)*4-1;
else
    first = 1992;
    last = 2005;
    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])';
f0data = reshape(f0data,[T,N])';
f1data = reshape(f1data,[T,N])';

clear Y
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));

rhohatvec(i,1) = xvec_ME;

%% Then, take 1992-2005 estimated parameters and use compute encompassing weights
dat2 = readmatrix(strcat('processed_data/spf9219_',macrovar(i),'_oos2.csv'));

id2 = dat2(:,1);
year2 = dat2(:,2);
quarter2 = dat2(:,3);
fe2 = dat2(:,4);
fr2 = dat2(:,5);
lfe2 = dat2(:,6);
f0data2 = dat2(:,7);
f1data2 = dat2(:,8);

if macrovar(i)=="cpi"
    first = 2007;
    last = 2019;
    T2 = (last-first+1)*4;
elseif macrovar(i)=="emp"
    first = 2012;
    last = 2019;
    T2 = (last-first+1)*4;
else
    first = 2006;
    last = 2019;
    T2 = (last-first+1)*4;
end

N2 = size(dat2(:,1),1)/T2;
fe2 = reshape(fe2,[T2,N2])';
fr2 = reshape(fr2,[T2,N2])';
lfe2 = reshape(lfe2,[T2,N2])';
f0data2 = reshape(f0data2,[T2,N2])';
f1data2 = reshape(f1data2,[T2,N2])';

Y2(:,:,1)=fe2;
Y2(:,:,2)=fr2;
Y2(:,:,3)=lfe2;


%% Report encompassing weights 

% run KF using parameter estimates
xvec_RE_all = [xvec(2:end) xvec_RE];
xvec_OC_all = [xvec(2:end) xvec_RE xvec_OC];
xvec_DE_all = [xvec(2:end) xvec_RE xvec_DE];
xvec_ME_all = [xvec(2:end) xvec_RE xvec_ME];

[Yact_RE,Yhat_RE,Yf1_RE] = predictfnc_RE(Y2,xvec_RE_all); 
[Yact_OC,Yhat_OC,Yf1_OC] = predictfnc_OC(Y2,xvec_OC_all); 
[Yact_DE,Yhat_DE,Yf1_DE] = predictfnc_DE(Y2,xvec_DE_all); 
[Yact_ME,Yhat_ME,Yf1_ME] = predictfnc_ME(Y2,xvec_ME_all); 

% Extract forecast (data vs. model)
f1_RE = Yf1_RE{1}(1:3:end);
f1_OC = Yf1_OC{1}(1:3:end);
f1_DE = Yf1_DE{1}(1:3:end);
f1_ME = Yf1_ME{1}(1:3:end);

for t=2:T2
    
    f1_RE = [f1_RE ; Yf1_RE{t}(1:3:end)];
    f1_OC = [f1_OC ; Yf1_OC{t}(1:3:end)];
    f1_DE = [f1_DE ; Yf1_DE{t}(1:3:end)];
    f1_ME = [f1_ME ; Yf1_ME{t}(1:3:end)];
end

f1tmp = f1data2(:,1);
f1tmp(isnan(f1tmp))=[];
f1vec = f1tmp;
for t=2:T2
    f1tmp = f1data2(:,t);
    f1tmp(isnan(f1tmp))=[];
    f1vec = [f1vec;f1tmp];
end


%% Regressing f0 and f1 on actual f0 and f1
RE_o = f1_RE;
OC_o = f1_OC;
DE_o = f1_DE;
ME_o = f1_ME;
Act_o = f1vec;

lb = [-Inf zeros(1,4)];
ub = [Inf ones(1,4)];
Aeq = [0 ones(1,4)];
beq =  1;
cons = ones(size(RE_o,1),1);
[coef,resnorm,residual,exitflag,output,lambda] = lsqlin([cons RE_o OC_o DE_o ME_o],Act_o,[],[],Aeq,beq,lb,ub);

Encompassing_Weights(i,:) = coef(2:end);
end %i
[macrovar' Encompassing_Weights]

% Figure C1
abbreviation = ["RGDP" "PGDP" "CPI" "EMP" "UNEMP" "RCONS" "RFED" "RSL" "RNRES" "RRES" "TB3M" "GB10Y" "HOUS" "IP"];

rho1vec   = fund_parameter_estimates(:,2);
rho2vec   = fund_parameter_estimates(:,3);
rhosum = rho1vec+rho2vec;
Encompassing_Weight_ME=Encompassing_Weights(:,4);
ME_bias = rhohatvec - rho1vec;
FigureC1_data = [abbreviation' macrovar' rhosum Encompassing_Weight_ME ME_bias];
writematrix(FigureC1_data,'/Users/julioortiz/Dropbox/ME/replication/results/FigureC1.xls')



