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

%Loop through subsamples
Nsamp=2;

xvec_RE = NaN(Nsamp,1);
xvec_OC = NaN(Nsamp,1);
xvec_DE = NaN(Nsamp,1);
xvec_ME = NaN(Nsamp,1);

fval_RE = NaN(Nsamp,1);
fval_OC = NaN(Nsamp,1);
fval_DE = NaN(Nsamp,1);
fval_ME = NaN(Nsamp,1);

AIC_RE = NaN(Nsamp,1);
AIC_OC = NaN(Nsamp,1);
AIC_DE = NaN(Nsamp,1);
AIC_ME = NaN(Nsamp,1);

subsample = 1:Nsamp;

%% First, solve actuals
dat=readmatrix('processed_data/rgdp_actuals.csv');

act_rt=dat(:,3);
tStart=tic;
[xvec_o,fval]=filterfnc_ar2_H(act_rt);
xvec_y=xvec_o;
tEnd=toc(tStart);
fprintf('Fundamental parameters estimated in %d minutes and %f seconds.\n',floor(tEnd/60), rem(tEnd,60));

%% Next, solve panel MLE 
daty=readmatrix('processed_data/spf9219_rgdp_3obs_subsample_young.csv');
dato=readmatrix('processed_data/spf9219_rgdp_3obs_subsample_old.csv');

Ty = (2019-1993+1)*4;
To = (2019-1992+1)*4;

Ny = size(daty(:,1),1)/Ty;
No = size(dato(:,1),1)/To;

fey = daty(:,4);
fry = daty(:,5);
lfey = daty(:,6);
f0datay = daty(:,7);
f1datay = daty(:,8);
fey=reshape(fey,[Ty,Ny])';
fry=reshape(fry,[Ty,Ny])';
lfey=reshape(lfey,[Ty,Ny])';
f0datay=reshape(f0datay,[Ty,Ny])';
f1datay=reshape(f1datay,[Ty,Ny])';

feo = dato(:,4);
fro= dato(:,5);
lfeo = dato(:,6);
f0datao = dato(:,7);
f1datao = dato(:,8);
feo=reshape(feo,[To,No])';
fro=reshape(fro,[To,No])';
lfeo=reshape(lfeo,[To,No])';
f0datao=reshape(f0datao,[To,No])';
f1datao=reshape(f1datao,[To,No])';

Yy(:,:,1)=fey;
Yy(:,:,2)=fry;
Yy(:,:,3)=lfey;

Yo(:,:,1)=feo;
Yo(:,:,2)=fro;
Yo(:,:,3)=lfeo;

tStart=tic;
[xvec_RE_y,fval_RE_y]=filtervecfnc_RE(Yy,Ny,Ty,xvec_y(2:4)); %RE model
[xvec_OC_y,fval_OC_y]=filtervecfnc_OC(Yy,Ny,Ty,xvec_y(2:4),xvec_RE_y); %OC model
[xvec_DE_y,fval_DE_y]=filtervecfnc_DE(Yy,Ny,Ty,xvec_y(2:4),xvec_RE_y); %DE model
[xvec_ME_y,fval_ME_y]=filtervecfnc_ME(Yy,Ny,Ty,xvec_y(2:4),xvec_RE_y); %ME model

[xvec_RE_o,fval_RE_o]=filtervecfnc_RE(Yo,No,To,xvec_o(2:4)); %RE model
[xvec_OC_o,fval_OC_o]=filtervecfnc_OC(Yo,No,To,xvec_o(2:4),xvec_RE_o); %OC model
[xvec_DE_o,fval_DE_o]=filtervecfnc_DE(Yo,No,To,xvec_o(2:4),xvec_RE_o); %DE model
[xvec_ME_o,fval_ME_o]=filtervecfnc_ME(Yo,No,To,xvec_o(2:4),xvec_RE_o); %ME model
tEnd=toc(tStart);
fprintf('Info/bias parameters estimated in %d minutes and %f seconds.\n',floor(tEnd/60), rem(tEnd,60));

%% Computing encompassing weights

% run KF using parameter estimates
xvec_RE_all_y = [xvec_y(2:end) xvec_RE_y];
xvec_OC_all_y = [xvec_y(2:end) xvec_RE_y xvec_OC_y];
xvec_DE_all_y = [xvec_y(2:end) xvec_RE_y xvec_DE_y];
xvec_ME_all_y = [xvec_y(2:end) xvec_RE_y xvec_ME_y];

[Yact_RE_y,Yhat_RE_y,Yf1_RE_y] = predictfnc_RE(Yy,xvec_RE_all_y); 
[Yact_OC_y,Yhat_OC_y,Yf1_OC_y] = predictfnc_OC(Yy,xvec_OC_all_y); 
[Yact_DE_y,Yhat_DE_y,Yf1_DE_y] = predictfnc_DE(Yy,xvec_DE_all_y); 
[Yact_ME_y,Yhat_ME_y,Yf1_ME_y] = predictfnc_ME(Yy,xvec_ME_all_y); 

xvec_RE_all_o = [xvec_o(2:end) xvec_RE_o];
xvec_OC_all_o = [xvec_o(2:end) xvec_RE_o xvec_OC_o];
xvec_DE_all_o = [xvec_o(2:end) xvec_RE_o xvec_DE_o];
xvec_ME_all_o = [xvec_o(2:end) xvec_RE_o xvec_ME_o];

[Yact_RE_o,Yhat_RE_o,Yf1_RE_o] = predictfnc_RE(Yo,xvec_RE_all_o); 
[Yact_OC_o,Yhat_OC_o,Yf1_OC_o] = predictfnc_OC(Yo,xvec_OC_all_o); 
[Yact_DE_o,Yhat_DE_o,Yf1_DE_o] = predictfnc_DE(Yo,xvec_DE_all_o); 
[Yact_ME_o,Yhat_ME_o,Yf1_ME_o] = predictfnc_ME(Yo,xvec_ME_all_o); 

% Extract forecasts (data and model, young)
f1_RE_y = Yf1_RE_y{1}(1:3:end);
f1_OC_y = Yf1_OC_y{1}(1:3:end);
f1_DE_y = Yf1_DE_y{1}(1:3:end);
f1_ME_y = Yf1_ME_y{1}(1:3:end);

for t=2:Ty
 
    f1_RE_y = [f1_RE_y ; Yf1_RE_y{t}(1:3:end)];
    f1_OC_y = [f1_OC_y ; Yf1_OC_y{t}(1:3:end)];
    f1_DE_y = [f1_DE_y ; Yf1_DE_y{t}(1:3:end)];
    f1_ME_y = [f1_ME_y ; Yf1_ME_y{t}(1:3:end)];
end

f1tmp_y = f1datay(:,1);
f1tmp_y(isnan(f1tmp_y))=[];
f1vec_y = f1tmp_y;
for t=2:Ty
    f1tmp_y = f1datay(:,t);
    f1tmp_y(isnan(f1tmp_y))=[];
    f1vec_y = [f1vec_y;f1tmp_y];
end

RE_o_y = f1_RE_y;
OC_o_y = f1_OC_y;
DE_o_y = f1_DE_y;
ME_o_y = f1_ME_y;
Act_o_y = f1vec_y;

lb    = [-Inf zeros(1,4)];
ub    = [Inf ones(1,4)];
Aeq   = [0 ones(1,4)];
beq   =  1;
cons  = ones(size(RE_o_y,1),1);
[coef_y,resnorm,residual,exitflag,output,lambda] = lsqlin([cons RE_o_y OC_o_y DE_o_y ME_o_y],Act_o_y,[],[],Aeq,beq,lb,ub);


% Extract forecasts (data and model, old)
f1_RE_o = Yf1_RE_o{1}(1:3:end);
f1_OC_o = Yf1_OC_o{1}(1:3:end);
f1_DE_o = Yf1_DE_o{1}(1:3:end);
f1_ME_o = Yf1_ME_o{1}(1:3:end);

for t=2:To
    f1_RE_o = [f1_RE_o ; Yf1_RE_o{t}(1:3:end)];
    f1_OC_o = [f1_OC_o ; Yf1_OC_o{t}(1:3:end)];
    f1_DE_o = [f1_DE_o ; Yf1_DE_o{t}(1:3:end)];
    f1_ME_o = [f1_ME_o ; Yf1_ME_o{t}(1:3:end)];
end

f1tmp_o = f1datao(:,1);
f1tmp_o(isnan(f1tmp_o))=[];
f1vec_o = f1tmp_o;
for t=2:To
    f1tmp_o = f1datao(:,t);
    f1tmp_o(isnan(f1tmp_o))=[];
    f1vec_o = [f1vec_o;f1tmp_o];
end

RE_o_o = f1_RE_o;
OC_o_o = f1_OC_o;
DE_o_o = f1_DE_o;
ME_o_o = f1_ME_o;
Act_o_o = f1vec_o;

lb    = [-Inf zeros(1,4)];
ub    = [Inf ones(1,4)];
Aeq   = [0 ones(1,4)];

beq   =  1;
cons  = ones(size(f1_RE_o,1),1);
[coef_o,resnorm,residual,exitflag,output,lambda] = lsqlin([cons RE_o_o OC_o_o DE_o_o ME_o_o],Act_o_o,[],[],Aeq,beq,lb,ub);

[coef_y';coef_o']
