%STDERR   compute bootstrapped standard errors.
%

%         Ellen McGrattan, 7-1-16
%         Revised, ERM, 1-15-17

load ../MLEdata/ldhpctot.dat
load ../MLEdata/ldgypc.dat
load ../MLEdata/ldhpc21.dat
load ../MLEparam/fixedpar.dat

global ZVAR
ZVAR  = [ldgypc(149:268,[1:15]),ldhpc21(149:268,[20,17,19])];

param=[ ...
   0.973987380000000
   0.983936860000000
   0.995000000000000
   0.987364760000000
   0.995000000000000
   0.976335760000000
   0.963338200000000
   0.957058380000000
   0.988589810000000
   0.138223740000000
   0.330325780000000
   0.325550860000000
   0.037926150000000
   0.071745700000000
   0.035687300000000
  -0.027167120000000
   0.026517240000000
   0.054167870000000
   0.038989980000000
   0.021903100000000
  -0.010898110000000
   0.027624960000000
   0.036794750000000
   0.080336570000000
   0.051403580000000
   0.014711900000000
  -0.010023020000000
   0.008438070000000
  -0.006948480000000
  -1.302693830000000
  -1.197412250000000
  -2.873862620000000
   1.427087270000000
   0.941599070000000
   0.712428840000000
   0.761068240000000
   1.167855590000000
   1.547388380000000
   2.201757850000000
   0.702501190000000
   0.191947680000000
   0.542150540000000
   0.611146300000000];

%
% Bootstrapped standard errors 
%

[L,Lt,ut,X0,Cbar,A,K,D] = mleseq(param,fixedpar);
B         = 500;   % number of bootstrap replications
T         = length(ut);
Ybar      = 0*ut;
nx        = length(X0);
Xt        = zeros(T+1,nx);
Data      = zeros(18*(T+1),B);
Y         = zeros(T+1,18);
Y(1,:)    = ZVAR(1,:);
Theta     = zeros(B,length(param));
thet0     = param;
rand('state',100462)
for i=1:B;
  disp(sprintf('Iteration = %g',i));
  %
  % Draw u's uniformly from sample {ut(1),ut(2)....ut(T)}
  % 
  Xt(1,:)   = X0';
  for j=1:T;
    k             = ceil(rand*T);  % draw number between 1 and T
    Ybar(j,:)     = Xt(j,:)*Cbar'+ut(k,:);
    Xt(j+1,:)     = Xt(j,:)*A'+ut(k,:)*K';
  end;
  for j=2:T+1;
    Y(j,:) = Y(j-1,:)*D'+Ybar(j-1,:);
  end;
  ZVAR      = Y;
  Data(:,i) = ZVAR(:);
  theta     = uncmin(thet0,'mleq',fixedpar);
  Theta(i,:)= theta';
  if (i>1);
    disp('  Estimates and deviations:')
    disp([thet0,std(Theta(1:i,:))']);
  end;
end;
se  = std(Theta)';

%
% Print results 
%

disp('Results from Maximum Likelihood Estimation')
disp('------------------------------------------')
disp(' ')
disp('[Theta, Standard Errors] ')
disp(sprintf(' %10.3e %10.3e\n', [param,se]'))
disp(' ')
fprintf('L(Theta) = %g ',L)
disp(' ')

