clear; clc; close all;

%set directory
cd('/Users/julioortiz/Dropbox/ME/replication/');
addpath '/Users/julioortiz/Dropbox/ME/replication/code/functions';

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

Encompassing_Weights = zeros(length(macrovar),4);

xmat    = zeros(length(macrovar),4);
xmat_se = zeros(length(macrovar),4);

for i=1:length(macrovar)
clear fe fr lfe Y 

%% First, solve for realized macro parameters
dat_actrt=readmatrix(strcat('processed_data/',macrovar(i),'_actuals.csv'));
act_rt=dat_actrt(:,3);
tStart=tic;
[xvec,fval,g,H]=filterfnc_ar2_H(act_rt);
xmat(i,:) = xvec;
tEnd=toc(tStart);
fprintf('Parameters estimated in %d minutes and %f seconds.\n',floor(tEnd/60), rem(tEnd,60));

%% Second, 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);
f0data=dat(:,8);
f1data=dat(:,9);

first = 1992;
last  = 2019;

if macrovar(i)=="cpi"
    first = 1994;
elseif macrovar(i)=="emp"
    first = 2004;
end

T = (last-first+1)*4;

if macrovar(i)=="tbond"
    T = (last-first+1)*4-1;
elseif macrovar(i)=="cpi"
    T = (last-first)*4+2;
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));

%% Third, calculate 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,Yf0_RE,Ycq_RE] = predictfnc_RE(Y,xvec_RE_all); 
[Yact_OC,Yhat_OC,Yf1_OC,Yf0_OC,Ycq_OC] = predictfnc_OC(Y,xvec_OC_all); 
[Yact_DE,Yhat_DE,Yf1_DE,Yf0_DE,Ycq_DE] = predictfnc_DE(Y,xvec_DE_all); 
[Yact_ME,Yhat_ME,Yf1_ME,Yf0_ME,Ycq_ME] = predictfnc_ME(Y,xvec_ME_all); 

% Extract model-implied data
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:T    
    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

% Extract actual data
f1data=reshape(f1data,[T,N])';
f1tmp = f1data(:,1);
f1tmp(isnan(f1tmp))=[];
f1vec = f1tmp;
for t=2:T
    f1tmp = f1data(:,t);
    f1tmp(isnan(f1tmp))=[];
    f1vec = [f1vec;f1tmp];
end %t

%Run constrained linear regression to obtain encompassing weights
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(f1_RE,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
Encompassing_Weights(i,:) = coef(2:end);

end %i


[macrovar' Encompassing_Weights]
