clear; clc; close all;

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

%% First, solve for realized macro parameters
dat_actrt=readmatrix('processed_data/rgdp_actuals.csv');
act_rt=dat_actrt(:,3);
tStart=tic;
[xvec,fval,g,H]=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('processed_data/spf9219_rgdp_spell4_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;

T = (last-first+1)*4;
N = size(dat(:,1),1)/T;
B=2000;

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

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

%% Bootstrap
rng(2501);
bfctmat=randi(N,B,N);

Yb=NaN(B,N,T,3);
for bct=1:B
    Yb(bct,:,:,:) = Y(bfctmat(bct,:),:,:);
end

for bct=1:B
    Ytmp=squeeze(Yb(bct,:,:,:));
    
    [xvec_RE,~]=filtervecfnc_RE(Ytmp,N,T,xvec(2:4));
    [xvec_OC,~]=filtervecfnc_OC(Ytmp,N,T,xvec(2:4),xvec_RE);
    [xvec_DE,~]=filtervecfnc_DE(Ytmp,N,T,xvec(2:4),xvec_RE);
    [xvec_ME,~]=filtervecfnc_ME(Ytmp,N,T,xvec(2:4),xvec_RE);
            
    folder = '/Users/julioortiz/Dropbox/ME/replication/results';
    
    filename = fullfile(folder, 'bootstrap_RE.txt');
    fid = fopen(filename, 'a');
    if fid == -1, error('Cannot open file: %s', filename); end
    fprintf(fid, '%f\r\n', xvec_RE);
    fclose(fid);
    
    filename = fullfile(folder, 'bootstrap_OC.txt');
    fid = fopen(filename, 'a');
    if fid == -1, error('Cannot open file: %s', filename); end
    fprintf(fid, '%f\r\n', xvec_OC);
    fclose(fid);
    
    filename = fullfile(folder, 'bootstrap_DE.txt');
    fid = fopen(filename, 'a');
    if fid == -1, error('Cannot open file: %s', filename); end
    fprintf(fid, '%f\r\n', xvec_DE);
    fclose(fid);
    
    filename = fullfile(folder, 'bootstrap_ME.txt');
    fid = fopen(filename, 'a');
    if fid == -1, error('Cannot open file: %s', filename); end
    fprintf(fid, '%f\r\n', xvec_ME);
    fclose(fid);
end %bct

%bootstrap SEs
thetamat_RE=readmatrix('/Users/julioortiz/Dropbox/ME/replication/results/bootstrap_RE.txt');
thetamat_OC=readmatrix('/Users/julioortiz/Dropbox/ME/replication/results/bootstrap_OC.txt');
thetamat_DE=readmatrix('/Users/julioortiz/Dropbox/ME/replication/results/bootstrap_DE.txt');
thetamat_ME=readmatrix('/Users/julioortiz/Dropbox/ME/replication/results/bootstrap_ME.txt');

%% Point estimate
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));

n_data=sum(~isnan(fe),'all');

% Model selection exercise
AIC_RE = 2*5 - 2*(-fval_RE);
AIC_OC = 2*6 - 2*(-fval_OC);
AIC_DE = 2*6 - 2*(-fval_DE);
AIC_ME = 2*6 - 2*(-fval_ME);

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

%% Store
fund_point_estimates = xvec;
fund_std_error       = sqrt(diag(inv(H)))';

bias_point_estimates = [xvec_RE xvec_OC xvec_DE xvec_ME];
bias_std_error       = [std(thetamat_RE) std(thetamat_OC) std(thetamat_DE) std(thetamat_ME)];

logL = [-fval_RE -fval_OC -fval_DE -fval_ME];

AIC = [AIC_RE AIC_OC AIC_DE AIC_ME];

Encompassing_Weight = coef(2:end)';

