%MONTE_UR  Take ndraw realizations of shocks from RBC model with
%          unit root in technology estimated by MLE with US data.
%          Use them to get ndraw realizations of productivity and 
%          hours.  Apply SVAR filter to see if true impulse responses
%          can be detected.
                                                                                
%          Ellen McGrattan, 3-17-04
%          Revised, ERM, 4-15-05


% Initialize random seeds
randn('seed',1004)
rand('seed',912)

% RBC model parameters
gn       = (1.01)^(1/4)-1;
gz       = (1.02)^(1/4)-1;
beta     = .98^(1/4);
delta    = 1-(1-.06)^(1/4);
psi      = 2.5;
sigma    = 1.000001;
theta    = .33;
param    = [gn;gz;beta;delta;psi;sigma;theta];
pq       = 1;    
if pq==1; pu;  end;  % load Sbar,P,Q for dataset PU
if pq==2; fr1; end;  % load Sbar,P,Q for dataset FR1
if pq==3; fr2; end;  % load Sbar,P,Q for dataset F-R2
if pq==4; cev; end;  % load Sbar,P,Q for dataset CEV
if pq==5; gr;  end;  % load Sbar,P,Q for dataset GR

% SVAR settings
alpha    = 0;
p        = 4;
nt       = 101;
nrep     = 1000;
ndraw    = 1500;
ldata    = 180;
Prd_ir   = zeros(nt,ndraw);
Hrs_ir   = zeros(nt,ndraw);
Prd_ciub = zeros(nt,ndraw);
Hrs_ciub = zeros(nt,ndraw);
Prd_cilb = zeros(nt,ndraw);
Hrs_cilb = zeros(nt,ndraw);
rho      = zeros(ndraw,1);
adf      = zeros(ndraw,1);
q        = zeros(ndraw,1);
rho      = zeros(ndraw,1);
LR       = zeros(ndraw,1);
AIC      = zeros(ndraw,12);
SIC      = zeros(ndraw,12);
for I=1:ndraw;
  mdata                = drawq_ur(param,Sbar,P,Q,ldata);
  [IR,CIlb,CIub]       = impreschk(mdata,alpha,p,nt,nrep);
  Prd_ir(:,I)          = IR(:,1);
  Hrs_ir(:,I)          = IR(:,2);
  Prd_ciub(:,I)        = CIub(:,1);
  Hrs_ciub(:,I)        = CIub(:,2);
  Prd_cilb(:,I)        = CIlb(:,1);
  Hrs_cilb(:,I)        = CIlb(:,2);
  [rho(I),adf(I),q(I)] = adftest(mdata(:,2),p,0,1);
  LR(I)                = lr(mdata,p,alpha);
  AIC(I,:)             = aic(mdata,12,alpha);
  SIC(I,:)             = sic(mdata,12,alpha);
end;
sHrsir   = sort(Hrs_ir');
sPrdir   = sort(Prd_ir');
mHrsir   = mean(sHrsir);
mPrdir   = mean(sPrdir);
mHrsciub = mean(Hrs_ciub');
mPrdciub = mean(Prd_ciub');
mHrscilb = mean(Hrs_cilb');
mPrdcilb = mean(Prd_cilb');

if ndraw>1
  lb   = round(.025*ndraw);
  ub   = ndraw-round(.025*ndraw);
  figure(1)
  plot([sPrdir(lb,:)',mPrdir',sPrdir(ub,:)'])
  title('Monte Carlo: Impulse reponse of productivity to technology shock')
  figure(2)
  plot([sHrsir(lb,:)',mHrsir',sHrsir(ub,:)'])
  title('Monte Carlo: Impulse reponse of hours to technology shock')
else
  figure(1)
  plot(Prd_ir)
  title('Impulse reponse of productivity to technology shock')
  figure(2)
  plot(Hrs_ir)
  title('Impulse reponse of hours to technology shock')
end;

