%Intang      Use parameters from intang_steady and first-order
% _Wedges    conditions from a 3-sector business cycle model 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 intang_steady.m first)
%

intang_params

bhat       = beta*(1+gamma)^(-sigma);
lshare1    = 1-thet1-phi1;
lshare2    = 1-thet2-phi2;
grate      = (1+eta)*(1+gamma);
maxit      = 100;

%
% 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;
xT         = xm-xn;

%
% Back out variables
%
g          = ym-xm-c;
kT         = kTs*ones(T,1);
for i=1:T-1;
  kTn(i)   = ((1-deltT)*kT(i)+xT(i))/grate;
  kT(i+1)  = kTn(i);
end;
kTn(T)     = ((1-deltT)*kT(T)+xT(T))/grate;
kTn        = kTn(:);
h1         = lshare1*(1-gtauh).*yb.*(1-h)./(psi*(1+tauc).*c);
h2         = h-h1-hn;
uc         = c.^(-sigma).*(1-h).^(psi*(1-sigma)); 
qxI        = lshare1/lshare2 * yb.*h2./h1;
k1T        = thet1*yb.*kT./(thet1*yb+thet2*qxI);
k2T        = kT-k1T;
xtau       = (1+gtaux).*(1-taud)./(1+tauc);
dtau       = (1-taud).*taup./(1+tauc);
ptau       = (1-taup).*(1-taud)./(1+tauc);
gtau       = (taug+(1-chi).*(1-gtauh)+chi.*(1-taup).*(1-taud))./(1+tauc);
ucn        = [uc(2:T);uc(T)];
ybn        = [yb(2:T);yb(T)];
qxIn       = [qxI(2:T);qxI(T)];
ptaun      = [ptau(2:T);ptau(T)];
gtaun      = [gtau(2:T);gtau(T)];
q          = ones(T,1);
kI         = kIs*ones(T,1);
kIn        = kI;
for it=1:maxit;

  xI       = qxI./q;
  for i=1:T-1;
    kIn(i) = ((1-deltI)*kI(i)+xI(i))/grate;
    kI(i+1)= kIn(i);
  end;
  kIn(T)   = ((1-deltI)*kI(T)+xI(T))/grate;
  qn       = [q(2:T);q(T)];
  res      = q.*gtau.*uc-bhat*ucn.* ...
                ((1-deltI)*qn.*gtaun+ptaun.*(phi1*ybn+phi2*qxIn)./kIn);
  dres     = zeros(T);

  for j=1:T;
    qp     = q;
    qp(j)  = q(j)+.0001;
 
    xI     = qxI./qp;
    for i=1:T-1;
      kIn(i)  = ((1-deltI)*kI(i)+xI(i))/grate;
      kI(i+1) = kIn(i);
    end;
    kIn(T) = ((1-deltI)*kI(T)+xI(T))/grate;
    qn     = [qp(2:T);qp(T)];
    tem    = qp.*gtau.*uc-bhat*ucn.* ...
                ((1-deltI)*qn.*gtaun+ptaun.*(phi1*ybn+phi2*qxIn)./kIn);
    dres(:,j) = (tem-res)/.0001;
  end;
  q0       = q;
  q        = q-dres\res;
end; 
if sqrt((q-q0)'*(q-q0))/T>1e-7;
  disp('Fixed point for q did not converge')
end;
xI         = qxI./q;
for i=1:T-1;
  kIn(i)   = ((1-deltI)*kI(i)+xI(i))/grate;
  kI(i+1)  = kIn(i);
end;
kIn(T)     = ((1-deltI)*kI(T)+xI(T))/grate;
z1         = (yb./(k1T.^thet1.*kI.^phi1)).^(1/lshare1)./h1;
z2         = (xI./(k2T.^thet2.*kI.^phi2)).^(1/lshare2)./h2;
a1         = z1.^lshare1;
a2         = z2.^lshare2;

%xwedge forwards
xwedge     = ones(T,1);
rT         = thet1*yb./k1T;
for i=1:T-1;

  rhs1        = bhat*uc(i+1)*(ptau(i+1)*(rT(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);

end;


inputs = [z1,z2,xwedge,g,tauc,gtauh,gtaux,tauk,taup,taud,taug,chi,yn,xn,hn];
initial= [kTn,kIn,h1,h2];

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