%RIntang     Use parameters from rintang_steady and first-order
% _Wedges    conditions from a 2-sector business cycle model with
%            intangible capital 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 rintang_steady.m first)
%

rintang_params

bhat       = beta*(1+gamma)^(-sigma);
lshare     = 1-theta-phi;
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);
taug       = taugs*ones(T,1);
chi        = chis*ones(T,1);
yb         = ym-yn;
xb         = xm-xn;
hb         = h-hn;

%
% Back out variables
%
g          = ym-xm-c;
xu         = psi*(1+tauc).*c.*hb./(lshare*(1-gtauh).*(1-h))-yb;
tyb        = yb+xI;
kT         = kTs*ones(T,1);
kI         = kIs*ones(T,1);
for i=1:T-1;
  kTn(i)   = ((1-deltT)*kT(i)+xT(i))/grate;
  kT(i+1)  = kTn(i);
  kIn(i)   = ((1-deltI)*kI(i)+xI(i))/grate;
  kI(i+1)  = kIn(i);
end;
kTn(T)     = ((1-deltT)*kT(T)+xT(T))/grate;
kTn        = kTn(:);
kIn(T)     = ((1-deltI)*kI(T)+xI(T))/grate;
kIn        = kIn(:);
uc         = c.^(-sigma).*(1-h).^(psi*(1-sigma)); 
xtau       = (1+gtaux).*(1-taud)./(1+tauc);
dtau       = (1-taud).*taup./(1+tauc);
ptau       = (1-taup).*(1-taud)./(1+tauc);
gtau       = ((1-chi).*(1-gtauh)+chi.*(1-taup).*(1-taud))./(1+tauc);
ucn        = [uc(2:T);uc(T)];
ybn        = [yb(2:T);yb(T)];
xun        = [xu(2:T);xu(T)];
tybn       = ybn+xun;
ptaun      = [ptau(2:T);ptau(T)];
gtaun      = [gtau(2:T);gtau(T)];

A          = tyb./(kT.^theta.*kI.^phi.*hb.^lshare);
xwedge     = ones(T,1);
gwedge     = ones(T,1);
rT         = theta*tyb./kT;
for i=1:T-1;

  rhs1        = bhat*uc(i+1)*(ptau(i+1)*(rb(i+1)-tauk(i+1))+deltT*dtau(i+1));
  rhs2        = bhat*(1-deltT)*uc(i+1)*xtau(i+1);
  xwedge(i+1) = rhs2/(uc(i)*xtau(i)/xwedge(i)-rhs1);
  gwedge(i+1) = (gwedge(i)*uc(i)+gtau(i)*uc(i)- ...
                   ((1-deltI)*bhat*uc(i+1)*gtau(i+1)+ ...
                    phi*bhat*uc(i+1)*ptau(i+1)*tyb(i+1)/kI(i+1)))/ ...
                  ((1-deltI)*bhat*uc(i+1));
end;
xwedge = xwedge(:);

inputs = [A,ones(T,1),xwedge,g,tauc,gtauh,gtaux,tauk,taup,taud, ...
          gwedge.*(1+tauc),chi,yn,xn,hn];
initial= [kTn,kIn,hb];

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

