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

%% First, estimate actuals
dat_actrt=readmatrix('rgdp_actuals.csv');
year=dat_actrt(:,1);
quarter=dat_actrt(:,2);
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));

%% Then estimate bias and info params
dat=readmatrix('processed_data/spf9219_rgdp_spell8_fefrlfe.csv');

B=2000;

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;

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

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

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=0.0;
    [xvec_DE,thetasol_DE]=filtervecfnc_DE2(Ytmp,N,T,xvec(2:4),xvec_RE);
    [xvec_ME,thetasol_ME]=filtervecfnc_ME2(Ytmp,N,T,xvec(2:4),xvec_RE);
            
    folder = '/Users/julioortiz/Dropbox/ME/replication/results';
        
    filename = fullfile(folder, 'bootstrap_DE2.txt');
    fid = fopen(filename, 'a');
    if fid == -1, error('Cannot open file: %s', filename); end
    fprintf(fid, '%f %f \n', xvec_DE);%[sigv varphi]
    fclose(fid);
    
    filename = fullfile(folder, 'bootstrap_ME2.txt');
    fid = fopen(filename, 'a');
    if fid == -1, error('Cannot open file: %s', filename); end
    fprintf(fid, '%f %f \n', xvec_ME);%[sigv rhohat]
    fclose(fid);
 end

%% Reporting full MLE estimates
%point estimate
[xvec_DE,thetasol_DE]=filtervecfnc_DE2(Y,N,T,xvec(2:4),xvec_RE);
[xvec_ME,thetasol_ME]=filtervecfnc_ME2(Y,N,T,xvec(2:4),xvec_RE);

%bootstrap SEs
thetamat_DE=readmatrix('/Users/julioortiz/Dropbox/ME/replication/results/bootstrap_DE2.txt');
thetamat_ME=readmatrix('/Users/julioortiz/Dropbox/ME/replication/results/bootstrap_ME2.txt');

%report results
[xvec_DE;std(thetamat_DE(:,1:end-1),[],1)]
[xvec_ME;std(thetamat_ME(:,1:end-1),[],1)]




