function [res,var,var2] = res_fullmodel(Z,param);
%RES_FULLMODEL  residual for multi-sector RBC model with intangible capital
%

%               Ellen R. McGrattan, 11-27-13
%               Revised, ERM, 11-27-13

%______________________________________________________________________________
%
% PARAMETERS


n          = param(1);
alpha      = param(2);
beta       = param(3);
deli       = param(5:4+n);
delt       = param(5+n:4+2*n);
gamma      = zeros(n);
gamma(:)   = param(5+3*n:4+n*n+3*n);
gn         = param(5+n*n+3*n);
gz         = param(6+n*n+3*n);
nu         = zeros(n);
nu(:)      = param(7+n*n+3*n:6+2*n*n+3*n);
omega      = param(7+2*n*n+3*n:6+2*n*n+4*n);
p1         = param(7+2*n*n+4*n);
phi        = param(8+2*n*n+4*n:7+2*n*n+5*n);
psi        = param(8+2*n*n+5*n);
rho        = param(9+2*n*n+5*n);
theta      = param(16+2*n*n+5*n:15+2*n*n+6*n);
zeta       = zeros(n);
zeta(:)    = param(16+2*n*n+8*n:15+3*n*n+8*n);
fkap       = param(17+3*n*n+8*n:16+3*n*n+9*n);
ftaub      = param(17+3*n*n+9*n:16+3*n*n+10*n);
nb         = param(18+3*n*n+11*n);
ib         = param(19+3*n*n+11*n:18+3*n*n+11*n+nb);

pvec       = param(19+3*n*n+11*n+nb:18+3*n*n+12*n+nb);
qvec       = param(19+3*n*n+12*n+nb:18+3*n*n+13*n+nb);
yvec       = param(19+3*n*n+13*n+nb:18+3*n*n+14*n+nb);
xvec       = param(19+3*n*n+14*n+nb:18+3*n*n+15*n+nb);
agghss     = param(19+3*n*n+15*n+nb);

h25ss      = param(20+3*n*n+15*n+nb);
h29ss      = param(21+3*n*n+15*n+nb);
h211ss     = param(22+3*n*n+15*n+nb);
wgt        = param(23+3*n*n+15*n+nb:22+3*n*n+16*n+nb);
dbar       = param(23+3*n*n+16*n+nb:22+3*n*n+17*n+nb);

o          = ones(n,1);
betg       = beta*(1+gz)^(-alpha);
irate      = 1/betg-1;
gam        = sum(gamma)';
lshare     = 1-theta-phi-gam;

m1mat      = zeros(n);
m2mat      = zeros(n);
m1matn     = zeros(n);
m2matn     = zeros(n);

%______________________________________________________________________________
%
% VARIABLES APPEARING IN RESIDUALS

kt         = exp(Z(1:n));
ki         = exp(Z(n+1:2*n));
b          = exp(Z(2*n+1:2*n+nb));

ktn        = exp(Z(2*n+nb+1:3*n+nb));
kin        = exp(Z(3*n+nb+1:4*n+nb));
bn         = exp(Z(4*n+nb+1:4*n+2*nb));

ktnn       = exp(Z(4*n+2*nb+1:5*n+2*nb));
kinn       = exp(Z(5*n+2*nb+1:6*n+2*nb));
bnn        = exp(Z(6*n+2*nb+1:6*n+3*nb));

irate      = Z(6*n+3*nb+1);
iraten     = Z(6*n+3*nb+2);

z1         = exp(Z(6*n+3*nb+3:7*n+3*nb+2));
z2         = exp(Z(7*n+3*nb+3:8*n+3*nb+2));
gov        = exp(Z(8*n+3*nb+3:9*n+3*nb+2));
tauc       = Z(9*n+3*nb+3);
taud       = Z(9*n+3*nb+4);
tauh       = Z(9*n+3*nb+5);
tauk       = Z(9*n+3*nb+6);
taup       = Z(9*n+3*nb+7);
taux       = Z(9*n+3*nb+8);
z          = exp(Z(9*n+3*nb+9));
fshk       = exp(Z(9*n+3*nb+10:9*n+4*nb+9));
fscom      = exp(Z(9*n+4*nb+10));

z1n        = exp(Z(9*n+4*nb+11:10*n+4*nb+10));
z2n        = exp(Z(10*n+4*nb+11:11*n+4*nb+10));
govn       = exp(Z(11*n+4*nb+11:12*n+4*nb+10));
taucn      = Z(12*n+4*nb+11);
taudn      = Z(12*n+4*nb+12);
tauhn      = Z(12*n+4*nb+13);
taukn      = Z(12*n+4*nb+14);
taupn      = Z(12*n+4*nb+15);
tauxn      = Z(12*n+4*nb+16);
zn         = exp(Z(12*n+4*nb+17));
fshkn      = exp(Z(12*n+4*nb+18:12*n+5*nb+17));
fscomn     = exp(Z(12*n+5*nb+18));

c          = exp(Z(12*n+5*nb+19:13*n+5*nb+18));
h1         = exp(Z(13*n+5*nb+19:14*n+5*nb+18));
h2         = exp(Z(14*n+5*nb+19:15*n+5*nb+18));
m1mat(:)   = exp(Z(15*n+5*nb+19:15*n+5*nb+n*n+18));
m2mat(:)   = exp(Z(15*n+5*nb+n*n+19:15*n+5*nb+2*n*n+18));

d          = Z(15*n+5*nb+2*n*n+19:15*n+6*nb+2*n*n+18);
wedge      = zeros(n,1);
wedge(ib)  = Z(15*n+6*nb+2*n*n+19:15*n+7*nb+2*n*n+18);

cn         = exp(Z(15*n+7*nb+2*n*n+19:16*n+7*nb+2*n*n+18));
h1n        = exp(Z(16*n+7*nb+2*n*n+19:17*n+7*nb+2*n*n+18));
h2n        = exp(Z(17*n+7*nb+2*n*n+19:18*n+7*nb+2*n*n+18));
m1matn(:)  = exp(Z(18*n+7*nb+2*n*n+19:18*n+7*nb+3*n*n+18));
m2matn(:)  = exp(Z(18*n+7*nb+3*n*n+19:18*n+7*nb+4*n*n+18));
dn         = Z(18*n+7*nb+4*n*n+19:18*n+8*nb+4*n*n+18);
wedgen     = zeros(n,1);
wedgen(ib) = Z(18*n+8*nb+4*n*n+19:18*n+9*nb+4*n*n+18);

mmat       = m1mat+m2mat;
mmatn      = m1matn+m2matn;
prodxt     = (1+gn)*(1+gz)*ktn-(1-delt).*kt;
prodxtn    = (1+gn)*(1+gz)*ktnn-(1-delt).*ktn;
prodxi     = (1+gn)*(1+gz)*kin-(1-deli).*ki;
prodxin    = (1+gn)*(1+gz)*kinn-(1-deli).*kin;
h          = h1+h2;
hn         = h1n+h2n;
p          = p1*(c/c(1)).^(-1/rho).*omega/omega(1);
pn         = p1*(cn/cn(1)).^(-1/rho).*omega/omega(1);
aggc       = (sum(omega.*c.^((rho-1)/rho)))^(rho/(rho-1));
aggcn      = (sum(omega.*cn.^((rho-1)/rho)))^(rho/(rho-1));
aggl       = 1-sum(h);
aggln      = 1-sum(hn);
aggp       = sum(omega.^rho.*p.^(1-rho))^(1/(1-rho));
aggpn      = sum(omega.^rho.*pn.^(1-rho))^(1/(1-rho));
uc         = aggc^(-alpha)*aggl^(psi*(1-alpha));
ucn        = aggcn^(-alpha)*aggln^(psi*(1-alpha));
aggw       = psi*(1+tauc)*aggp*aggc/((1-tauh)*aggl);
aggwn      = psi*(1+taucn)*aggpn*aggcn/((1-tauhn)*aggln);
w          = aggw*o;
wn         = aggwn*o;
m1         = (prod(m1mat.^gamma)').^(o./gam);
m1n        = (prod(m1matn.^gamma)').^(o./gam);
m2         = (prod(m2mat.^gamma)').^(o./gam);
m2n        = (prod(m2matn.^gamma)').^(o./gam);
qxipy      = h2./h1;
qxipyn     = h2n./h1n;
k1t        = kt./(o+qxipy);
k1tn       = ktn./(o+qxipyn);
k2t        = kt-k1t;
k2tn       = ktn-k1tn;
zz         = (z*ones(n,1)).^wgt;
zz(n)      = 1;
zzn        = (zn*ones(n,1)).^wgt;
zzn(n)     = 1;
y          = k1t.^theta.*ki.^phi.*m1.^gam.*(zz.*z1.*h1).^lshare;
yn         = k1tn.^theta.*kin.^phi.*m1n.^gam.*(zzn.*z1n.*h1n).^lshare;
xif        = k2t.^theta.*ki.^phi.*m2.^gam.*(zz.*z2.*h2).^lshare;
xifn       = k2tn.^theta.*kin.^phi.*m2n.^gam.*(zzn.*z2n.*h2n).^lshare;
qxip       = qxipy.*y;
qxipn      = qxipyn.*yn;
qxi        = qxip.*p;
qxin       = qxipn.*pn;
xiy        = qxipy.^(1-phi).*(z2./z1).^lshare;
xiyn       = qxipyn.^(1-phi).*(z2n./z1n).^lshare;
q          = p.*qxipy./xiy;
qn         = pn.*qxipyn./xiyn;
pt         = prod(((p*o')./(zeta+1e-15)).^zeta)';
ptn        = prod(((pn*o')./(zeta+1e-15)).^zeta)';
pi         = prod(((q*o')./(nu+1e-15)).^nu)';
pin        = prod(((qn*o')./(nu+1e-15)).^nu)';
rtn        = theta.*pn.*yn./k1tn;
rin        = phi.*(pn.*yn+qxin)./kin;
xtmat      = zeta.*((o./p)*pt').*(o*prodxt');
xtmatn     = zeta.*((o./pn)*ptn').*(o*prodxtn');
ximat      = nu.*((o./q)*pi').*(o*prodxi');
ximatn     = nu.*((o./qn)*pin').*(o*prodxin');
m1matest   = gamma.*((o./p)*p').*(o*y').*(o*(1-wedge/(1-taup))');
m2matest   = gamma.*((o./p)*p').*(o*qxip').*(o*(1-wedge/(1-taup))');
m1matnest  = gamma.*((o./pn)*pn').*(o*yn').*(o*(1-wedgen/(1-taupn))');
m2matnest  = gamma.*((o./pn)*pn').*(o*qxipn').*(o*(1-wedgen/(1-taupn))');
lambda     = (1-taud)*uc/(aggp*(1+tauc))*o;
lambdan    = (1-taudn)*ucn/(aggpn*(1+taucn))*o;

lambda(ib) = lambda(ib)./(1+2*fkap(ib).*(d-dbar(ib)));
lambdan(ib)= lambdan(ib)./(1+2*fkap(ib).*(dn-dbar(ib)));
R          = 1+irate*(1-ftaub(ib));
Rn         = 1+iraten*(1-ftaub(ib));
defd       = (1-taup)*(p.*y+qxi-w.*h-sum((p*o').*mmat)'- ...
                sum((q*o').*ximat)')-(1+taux)*sum((p*o').*xtmat)'- ...
                (1-taup)*tauk*p.*kt+taup*delt.*p.*kt;
defd       = defd(ib)+bn./R-b-fkap(ib).*(d-dbar(ib)).^2;
defdn      = (1-taupn)*(pn.*yn+qxin-wn.*hn-sum((pn*o').*mmatn)'- ...
                sum((qn*o').*ximatn)')-(1+tauxn)*sum((pn*o').*xtmatn)'- ...
                (1-taupn)*taukn*pn.*ktn+taupn*delt.*pn.*ktn;
defdn      = defdn(ib)++bnn./Rn-bn-fkap(ib).*(dn-dbar(ib)).^2;
fshka      = zeros(n,1);
fshka(ib)  = fshk;


res        = [y-c-sum(xtmat')'-sum(m1mat')'-sum(m2mat')'-gov;
              yn-cn-sum(xtmatn')'-sum(m1matn')'-sum(m2matn')'-govn;
              xif-sum(ximat')';
              xifn-sum(ximatn')';
              aggw*o-(1-wedge/(1-taup)).*lshare.*p.*y./h1;
              aggwn*o-(1-wedgen/(1-taupn)).*lshare.*pn.*yn./h1n;
              m1mat(:)-m1matest(:);
              m1matn(:)-m1matnest(:);
              m2mat(:)-m2matest(:);
              m2matn(:)-m2matnest(:);
              fshk.*(p(ib).*ktn(ib)-bn/(1+irate))-p(ib).*y(ib)-q(ib).*xif(ib);
              fshkn.*(pn(ib).*ktnn(ib)-bnn/(1+iraten))-pn(ib).*yn(ib)-qn(ib).*xifn(ib);
              defd-d;
              defdn-dn;
              o-betg*lambdan.*pn.* ...
                    (1+(1-taupn)*(rtn./pn-delt-taukn) ...
                    -wedgen.*rtn./pn+((1+tauxn)*ptn./pn-1).*(1-delt))./ ...
                    (lambda.*((1+taux)*pt-wedge.*fshka.*p));
              o-betg*lambdan.* ...
                   ((1-taupn-wedgen).*rin+(1-taupn)*pin.*(1-deli))./ ...
                    ((1-taup)*lambda.*pi);
              o(ib)-R.*(betg*lambdan(ib)./lambda(ib)+wedge(ib).*fshka(ib)/(1+irate));
              1-betg*(1+irate)*ucn*aggp/(uc*aggpn)];


var        = [log((pvec(1)*y(1)+qvec(1)*xif(1))/(pvec(1)*yvec(1)+qvec(1)*xvec(1)));
              log((pvec(2)*y(2)+qvec(2)*xif(2))/(pvec(2)*yvec(2)+qvec(2)*xvec(2)));
              log((pvec(3)*y(3)+qvec(3)*xif(3))/(pvec(3)*yvec(3)+qvec(3)*xvec(3)));
              log((pvec(4)*y(4)+qvec(4)*xif(4))/(pvec(4)*yvec(4)+qvec(4)*xvec(4)));
              log(y(5)/yvec(5));
              log((pvec(6)*y(6)+qvec(6)*xif(6))/(pvec(6)*yvec(6)+qvec(6)*xvec(6)));
              log((pvec(7)*y(7)+qvec(7)*xif(7))/(pvec(7)*yvec(7)+qvec(7)*xvec(7)));
              log((pvec(8)*y(8)+qvec(8)*xif(8))/(pvec(8)*yvec(8)+qvec(8)*xvec(8)));
              log(y(9)/yvec(9));
              log((pvec(10)*y(10)+qvec(10)*xif(10))/(pvec(10)*yvec(10)+qvec(10)*xvec(10)));
              log(y(11)/yvec(11));
              log((pvec(12)*y(12)+qvec(12)*xif(12))/(pvec(12)*yvec(12)+qvec(12)*xvec(12)));
              log((pvec(13)*y(13)+qvec(13)*xif(13))/(pvec(13)*yvec(13)+qvec(13)*xvec(13)));
              log((pvec(14)*y(14)+qvec(14)*xif(14))/(pvec(14)*yvec(14)+qvec(14)*xvec(14)));
              log((pvec(15)*y(15)+qvec(15)*xif(15))/(pvec(15)*yvec(15)+qvec(15)*xvec(15)));
              log(h2(5)/h25ss);
              log(h2(9)/h29ss);
              log(h2(11)/h211ss);
];


%
% Additional variables:
%
tem = [ ...
   1.6373e-03   2.4804e-07
   2.4973e-03   9.6183e-06
   3.2746e-03   1.4075e-05
   1.2847e-02   2.5555e-05
   3.2775e-02   4.5600e-03
   1.8983e-02   5.6099e-04
   2.0330e-02   1.2720e-04
   1.0223e-02   1.9007e-05
   6.4827e-03   1.9918e-03
   2.9632e-02   5.2791e-04
   2.9987e-02   2.0725e-02
   3.7586e-02   3.5376e-04
   1.3937e-02   8.1012e-06
   1.1117e-02   1.3683e-04
   1.9256e-03   1.3647e-05];

h1ss = tem(:,1);
h2ss = tem(:,2);
var2 =[log((1-aggl)/agghss);
       log(sum(pvec.*y+qvec.*xif)/sum(pvec.*yvec+qvec.*xvec));
       log(pvec.*y+qvec.*xif);
       log((h1+h2)./(h1ss+h2ss));
       log(h1./h1ss);
       log(h2([5,9,11])./h2ss([5,9,11]))
       log(sum((pvec*ones(1,15)).*xtmat)');
       log(sum((qvec*ones(1,15)).*ximat)');
       log(pvec.*(sum(xtmat')'));
       log(qvec.*xif);
       sum(pvec.*(sum(xtmat')'))/sum(pvec.*yvec+qvec.*xvec);
       sum(qvec.*xif)/sum(pvec.*yvec+qvec.*xvec);
       wedge/(1-taup)];

%    =[log(total hours rel to trend)
%      log(total gross output rel to trend)
%      log(sectoral gross outputs)
%      log(sectoral hours)
%      log(hours in tangible production)
%      log(hours in intangible production)
%      log(tangible investment by using industry)
%      log(intangible investment expenditures by using industry)
%      log(tangible investment expenditures by commodity)
%      log(intangible investment expenditures by commodity)
%      log(tangible investment expenditures expenditures)
%      log(intangible investment expenditures)
%      implied labor wedges

