%------------------------------------------------------------------------------
% Bootstrap_2s Compute bootstrap standard errors for the 2-shock model
%              in mleq_ur_2s.m.
                                                                                
%              Ellen McGrattan, 2-16-02
%              Revised, ERM, 4-10-05
                                                                                
%------------------------------------------------------------------------------
load uszvarq_pu.dat
uszvarq_ur=uszvarq_pu;
global ZVAR
t    = uszvarq_ur(:,1);
ZVAR = uszvarq_ur(:,2:5);

[x1,f,g,code,status]   = uncmin(x0,'mleq_ur_2s',0);
[x1,f,g,code,status]   = uncmin(x1,'mleq_ur_2s',0);
[L,Sbar,P0,P,Q]        = mleq_ur_2s(x1);

ind   = [1,1,0,0, ...
         0,0,0,0, 0,1,0,0, 0,0,0,0, 0,0,0,0, ...
         1,0,0,0,   1,0,0,     0,0,       0]';
param = zeros(30,1);
param(find(ind)) = x1;

[L,Lt] = mleseq_ur(param,ind);

del = diag(max(abs(param)*1e-4,1e-8));
for i=1:length(param);
  if ind(i);
    [f1,f2] = mleseq_ur(param+del(:,i),ind);
    [m1,m2] = mleseq_ur(param-del(:,i),ind);
    dL(i,1)  = (f1-m1)/(2*del(i,i));
    dLt(i,:) = (f2-m2)'/(2*del(i,i));
  end;
end;
dLt =dLt(find(ind),:);
sum1=0;
[n,m]=size(dLt);
for t=1:m;
  sum1=sum1+dLt(:,t)*dLt(:,t)';
end;
se1 = diag(sqrt(inv(sum1)));

for i=1:length(param);
  if ind(i);
    for j=1:i;
      d2L(i,j)=(mleseq_ur(param+del(:,i)+del(:,j),ind)- ...
                mleseq_ur(param+del(:,i)-del(:,j),ind)- ...
                mleseq_ur(param-del(:,i)+del(:,j),ind)+ ...
                mleseq_ur(param-del(:,i)-del(:,j),ind))/(4*del(i,i)*del(j,j));
      if (i~=j); d2L(j,i)=d2L(i,j); end;
    end;
  end;
end;
d2L = d2L(find(ind),find(ind));
se2 = diag(sqrt(inv(d2L)));

[L,Lt,ut,X0,Cbar,A,K,D,gz] = mleseq_ur(param,0);
B         = 500;   % number of bootstrap replications
T         = length(ut);
Ybar      = 0*ut;
nx        = length(X0);
Xt        = zeros(T+1,nx);
Y         = zeros(T+1,4);
Y(1,:)    = log(ZVAR(1,:));
Data      = zeros(4*(T+1),B);
thet0     = param(find(ind));
Theta     = zeros(B,length(thet0));
rand('state',100462)
for i=1:B;
  %
  % 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      = exp(Y);
  Data(:,i) = ZVAR(:);
  theta     = uncmin(thet0,'mleq_ur_2s',0);
  theta(4:5)= abs(theta(4:5));
  Theta(i,:) = theta';
end;
se3 = std(Theta)';

disp('Results from Maximum Likelihood Estimation')
disp('------------------------------------------')
disp(' ')
disp('                          Standard Errors       ')
disp('                 -------------------------------')
disp('   Parameters      dL*dL''      d2L     bootstrap')
disp('   ---------------------------------------------')
disp(sprintf('   %9.3e     %9.3e  %9.3e  %9.3e\n',[x1,se1,se2,se3]'))
fprintf('   Log-Likelihood = %9.3f ',L)
disp(' ')
                                                                                

