%NoIntang    Use parameters from nointang_steady and first-order 
% _Wedges    conditions from a prototype 2-sector business cycle 
%            model (with a business and nonbusiness sector) to 
%            back out series 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-13-07

%---------------------------------------------------------------------
% Parameters (for steady state parameters, run nointang_steady.m first)
%

nointang_params

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

%
% Load BEA/FOF/CPS data 
%
load ../../DATA/observables.dat
t          = observables(:,1);
gdp        = observables(:,2)*gdps;
ym         = gdp;
T          = length(gdp);
yn         = observables(:,3)*gdps;
c          = observables(:,4)*gdps;
xm         = observables(:,5)*gdps;
xn         = observables(:,6)*gdps;
h          = observables(:,7);
hn         = observables(:,8);
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);
yb         = ym-yn;
xb         = xm-xn;
hb         = h-hn;

%
% Back out variables
%
g          = ym-xm-c;
kb         = kbs*ones(T,1);
for i=1:T-1;
  kbn(i)   = (delt1*kb(i)+xb(i))/grate;
  kb(i+1)  = kbn(i);
end;
kbn(T)     = (delt1*kb(T)+xb(T))/grate;
kbn        = kbn(:);
lwedge     = (psi*(1+tauc).*c.*hb)./(thet1*yb.*(1-h).*(1-gtauh));
A          = yb./(kb.^theta.*hb.^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)); 
rb         = theta*yb./kb;
%xwedge forwards
for i=1:T-1;

  rhs1        = bhat*uc(i+1)*(ptau(i+1)*(rb(i+1)-tauk(i+1))+delta*dtau(i+1));
  rhs2        = bhat*delt1*uc(i+1)*xtau(i+1);
  xwedge(i+1) = rhs2/(uc(i)*xtau(i)/xwedge(i)-rhs1);

end;

inputs     = [A,lwedge,xwedge,g,tauc,gtauh,gtaux,tauk,taup,taud,yn,xn,hn];
initial    = [kbn,h];
                                                                                
save inputs.dat inputs -ascii -double
save vinitial.dat initial -ascii -double
