%BCA_Wedges  Use parameters from bca_steady and first-order conditions
%            from a prototype business cycle model to back out the  
%            series which are used to form the wedges for business 
%            cycle accounting.  See details in the Appendix to:
%
%               Unmeasured Investment and the 1990s US Hours Boom
%               by Ellen McGrattan and Ed Prescott
                                                                                
%            Ellen McGrattan, 11-1-05
%            Revised, ERM, 8-6-07
                                                                                
%---------------------------------------------------------------------
% Parameters (for steady state parameters, run bca_steady.m first)
%

bca_params

bhat       = beta*(1+gamma)^(-sigma);
thet1      = 1-theta;
delt1      = 1-delta;
grate      = (1+eta)*(1+gamma);

% Load NIPA data 
%
load ../../DATA/observables.dat
t          = observables(:,1);
gdp        = observables(:,2)*gdps;
T          = length(gdp);
y          = gdp;
c          = observables(:,4)*gdps;
x          = observables(:,5)*gdps;
h          = observables(:,7);
tauc       = observables(:,9);
gtauh      = observables(:,11);
gtaux      = tauxs*ones(T,1);
tauk       = tauks*ones(T,1);
taup       = taups*ones(T,1);
taud       = tauds*ones(T,1);

%
% Back out variables
%
g          = gdp-x-c;
k          = ks*ones(T,1);
kn         = k;
for i=1:T-1;
  kn(i)    = (delt1*k(i)+x(i))/grate;
  k(i+1)   = kn(i);
end;
kn(T)      = (delt1*k(T)+x(T))/grate;
lwedge     = psi*(1+tauc).*c.*h./(thet1*y.*(1-h).*(1-gtauh));    
A          = y./(k.^theta.*h.^thet1);
xwedge     = ones(T,1);
xtau       = (1+gtaux).*(1-taud)./(1+tauc);
dtau       = (1-taud).*taup./(1+tauc);
ptau       = (1-taup).*(1-taud)./(1+tauc);
uc         = c.^(-sigma).*(1-h).^(psi*(1-sigma)); 
r          = theta*y./k;
%xwedge forwards
for i=1:T-1;
  xwedge(i+1) = (bhat*delt1*uc(i+1)*xtau(i+1))/ ...
                (uc(i)*xtau(i)/xwedge(i)- ...
                   bhat*uc(i+1)*ptau(i+1)*(r(i+1)-tauk(i+1))- ...
                   bhat*uc(i+1)*delta*dtau(i+1));

end;
%xwedge backards
%xwedge(T)  = uc(T)*xtau(T);
%for i=T-1:-1:1;
%  xwedge(i)  = bhat*delt1*xwedge(i+1)+ ...
%                 bhat*uc(i+1)*ptau(i+1)*(r(i+1)-tauk(i+1))+ ...
%                 bhat*uc(i+1)*delta*dtau(i+1);
%end;
%xwedge    = uc.*xtau./xwedge(:);

inputs     = [A,lwedge,xwedge,g,tauc,gtauh,gtaux,tauk,taup,taud];
initial    = [kn,h];

save inputs.dat inputs  -ascii -double
save vinitial.dat initial  -ascii -double

