rho      = -2;
sigma    = 1.5;
mtau2    = .0;
mtaun    = .25;
delu1    = .0;
Mxu1     = .0285;   

% From data 

Dn       = .25;                % assume working 25% of time
Dgamma   = .02/4;                % qtrly growth in technology
Deta     = .01/4;                % qtrly growth in population
Dc       = .5394;              % .576 - .087 (sub+ctx) +.0759 (kcd+kg ser)
Dxm1     = .0990;              % corporate fixed investment+.6*inventories
Dx2      = .1845;              % priv. inv.+gov. inv.+NE+durables-Dxm1
Dg       = .1901;              % government consumption 

Dw1      = .3780;              % corporate employee compensation
Dw2      = .2554;              % noncorporate employee compensation
Dr1x     = .1620;              % VA corp-W corp-.015-IBT corp
Dr2      = .2176;              % VA noncorp-W noncorp-.015+.076-IBT non 

Dkm1     = .95 *4;                % corporate tangible assets
Dk2      = 2.217 *4;              % tangible assets of HH, noncorp. + gov

Dibt     = .047;               % IBT 
mtau1    = .026/.073;          % corporate profit tax rate
mtauk1   = .02/Dkm1;           % property tax rate in corporate sector
mtauk2   = .012/Dk2;           % property tax rate in noncorporate sector


% Model parameters:

gamma    = Dgamma;
eta      = Deta;
mtauc    = Dibt/Dc;
Dku1     = Mxu1/(delu1+gamma+eta+gamma*eta);
delm1    = Dxm1/Dkm1 - (gamma+eta+gamma*eta);
del2     = Dx2/Dk2 - (gamma+eta+gamma*eta);
mphi     = Dr2/(Dw2+Dr2);
Mr2      = Dr2/Dk2;
beta     = (1+gamma)^sigma/((1-mtau2)*(Mr2-del2-mtauk2)+1);
betg     = beta*(1+gamma)^(-sigma);
Mru1     = (1-mtau2)*(Mr2-del2-mtauk2)+delu1;
Mrm1     = (Mru1-delu1)/(1-mtau1)+delm1+mtauk1;
Dy       = Dc+Dxm1+Mxu1+Dx2+Dg;
Dr1      = Dr1x+Mxu1;
Dva1     = Dr1+Dw1;
Dva2     = Dr2+Dw2;
mthetm   = Mrm1*Dkm1/Dva1;
mthetu   = Mru1*Dku1/Dva1;
Mw       = ((1-mphi)*Dva2+(1-mthetm-mthetu)*Dva1)/Dn;
Mn1      = Dw1/Mw;
Mn2      = Dw2/Mw;
psi      = (1-mtaun)*Mw*(1-Mn1-Mn2)/(1+mtauc)/Dc;
My1      = Dkm1^mthetm*Dku1^mthetu*Mn1^(1-mthetm-mthetu);
My2      = Dk2^mphi*Mn2^(1-mphi);
Mp1      = Dva1/My1;
Mp2      = Dva2/My2;
Aa1      = Mp1*(Dy/My1)^(rho-1);
Aa2      = Mp2*(Dy/My2)^(rho-1);
scale    = (Aa1+Aa2)^(1/rho);
a1       = Aa1/scale^rho;
a2       = Aa2/scale^rho;
grate    = gamma+eta+eta*gamma;
mg       = Dg;


%______________________________________________________________________________
%
% STEADY STATE FOR THE ECONOMY
%

y1y    = My1/Dy;  % Guess of corporate GNP to total GNP
n1     = Mn1;     % Guess of hours of work in the corporate sector

z0     = [y1y;n1];
del    = diag(z0) * .00001;

for i=1:20;

   y1y = z0(1); 
   n1  = z0(2); 
   y2y = ((scale^(-rho)-a1*y1y^rho)/a2)^(1/rho);
   p1  = scale^rho*a1*y1y^(rho-1);
   p2  = scale^rho*a2*y2y^(rho-1);
   r2  = (1/betg-1)/(1-mtau2)+del2+mtauk2;
   ru1 = 1/betg-1+delu1;
   rm1 = (1/betg-1)/(1-mtau1)+delm1+mtauk1;
   nkm1= (rm1/mthetm/p1 * (mthetm*ru1/mthetu/rm1)^mthetu )^ ...
         (1/(1-mthetm-mthetu)); 
   nk2 = (r2/mphi/p2)^(1/(1-mphi));
   km1 = n1/nkm1;
   ku1 = mthetu*rm1*km1/mthetm/ru1;
   nku1= n1/ku1;
   y1  = km1^mthetm*ku1^mthetu*n1^(1-mthetm-mthetu);
   y   = y1/y1y;
   y2  = (y^rho/(scale^rho*a2)-a1/a2*y1^rho)^(1/rho);
   k2  = y2*nk2^(mphi-1);
   n2  = nk2*k2; 
   xm1 = (grate+delm1)*km1;
   xu1 = (grate+delu1)*ku1;
   x2  = (grate+del2)*k2;
   c   = y-xm1-xu1-x2-mg;

   f   = [p1*(1-mthetm-mthetu)*nkm1^(-mthetm)*nku1^(-mthetu)- ...
             p2*(1-mphi)*nk2^(-mphi);
          psi*c-(1-mtaun)*p2*(1-mphi)*nk2^(-mphi)*(1-n1-n2)/(1+mtauc)];

   y1y = z0(1) + del(1,1);
   y2y = ((scale^(-rho)-a1*y1y^rho)/a2)^(1/rho);
   p1  = scale^rho*a1*y1y^(rho-1);
   p2  = scale^rho*a2*y2y^(rho-1);
   r2  = (1/betg-1)/(1-mtau2)+del2+mtauk2;
   ru1 = 1/betg-1+delu1;
   rm1 = (1/betg-1)/(1-mtau1)+delm1+mtauk1;
   nkm1= (rm1/mthetm/p1 * (mthetm*ru1/mthetu/rm1)^mthetu )^ ...
         (1/(1-mthetm-mthetu)); 
   nk2 = (r2/mphi/p2)^(1/(1-mphi));
   km1 = n1/nkm1;
   ku1 = mthetu*rm1*km1/mthetm/ru1;
   nku1= n1/ku1;
   y1  = km1^mthetm*ku1^mthetu*n1^(1-mthetm-mthetu);
   y   = y1/y1y;
   y2  = (y^rho/(scale^rho*a2)-a1/a2*y1^rho)^(1/rho);
   k2  = y2*nk2^(mphi-1);
   n2  = nk2*k2;
   xm1 = (grate+delm1)*km1;
   xu1 = (grate+delu1)*ku1;
   x2  = (grate+del2)*k2;
   c   = y-xm1-xu1-x2-mg;

   df(:,1) ...
       = ([p1*(1-mthetm-mthetu)*nkm1^(-mthetm)*nku1^(-mthetu)- ...
             p2*(1-mphi)*nk2^(-mphi);
           psi*c-(1-mtaun)*p2*(1-mphi)*nk2^(-mphi)*(1-n1-n2)/(1+mtauc)]- ...
           f)/del(1,1);

   y1y = z0(1);
   n1  = z0(2) + del(2,2);
   y2y = ((scale^(-rho)-a1*y1y^rho)/a2)^(1/rho);
   p1  = scale^rho*a1*y1y^(rho-1);
   p2  = scale^rho*a2*y2y^(rho-1);
   r2  = (1/betg-1)/(1-mtau2)+del2+mtauk2;
   ru1 = 1/betg-1+delu1;
   rm1 = (1/betg-1)/(1-mtau1)+delm1+mtauk1;
   nkm1= (rm1/mthetm/p1 * (mthetm*ru1/mthetu/rm1)^mthetu )^ ...
         (1/(1-mthetm-mthetu)); 
   nk2 = (r2/mphi/p2)^(1/(1-mphi));
   km1 = n1/nkm1;
   ku1 = mthetu*rm1*km1/mthetm/ru1;
   nku1= n1/ku1;
   y1  = km1^mthetm*ku1^mthetu*n1^(1-mthetm-mthetu);
   y   = y1/y1y;
   y2  = (y^rho/(scale^rho*a2)-a1/a2*y1^rho)^(1/rho);
   k2  = y2*nk2^(mphi-1);
   n2  = nk2*k2;
   xm1 = (grate+delm1)*km1;
   xu1 = (grate+delu1)*ku1;
   x2  = (grate+del2)*k2;
   c   = y-xm1-xu1-x2-mg;

   df(:,2) ...
       = ([p1*(1-mthetm-mthetu)*nkm1^(-mthetm)*nku1^(-mthetu)- ...
             p2*(1-mphi)*nk2^(-mphi);
           psi*c-(1-mtaun)*p2*(1-mphi)*nk2^(-mphi)*(1-n1-n2)/(1+mtauc)]- ...
           f)/del(2,2);

   z0  = z0 - df\f;

end;
y1y  = z0(1);
n1   = z0(2);
y2y  = ((scale^(-rho)-a1*y1y^rho)/a2)^(1/rho);
p1   = scale^rho*a1*y1y^(rho-1);
p2   = scale^rho*a2*y2y^(rho-1);
r2   = (1/betg-1)/(1-mtau2)+del2+mtauk2;
ru1  = 1/betg-1+delu1;
rm1  = (1/betg-1)/(1-mtau1)+delm1+mtauk1;
nkm1 = (rm1/mthetm/p1 * (mthetm*ru1/mthetu/rm1)^mthetu )^ ...
       (1/(1-mthetm-mthetu)); 
nk2  = (r2/mphi/p2)^(1/(1-mphi));
km1  = n1/nkm1;
ku1  = mthetu*rm1*km1/mthetm/ru1;
nku1 = n1/ku1;
k1   = km1+ku1;
y1   = km1^mthetm*ku1^mthetu*n1^(1-mthetm-mthetu);
y    = y1/y1y;
y2   = (y^rho/(scale^rho*a2)-a1/a2*y1^rho)^(1/rho);
k2   = y2*nk2^(mphi-1);
n2   = nk2*k2;
xm1  = (grate+delm1)*km1;
xu1  = (grate+delu1)*ku1;
x2   = (grate+del2)*k2;
x1   = xm1+xu1;
x    = x1+x2;
km1f = (1+gamma)*km1;
ku1f = (1+gamma)*ku1;
k1f  = km1f+ku1f;
k2f  = (1+gamma)*k2;
c    = y-xm1-xu1-x2-mg;
rf   = (1+gamma)^sigma/beta-1;
rm1  = p1*mthetm*y1/km1;
ru1  = p1*mthetu*y1/ku1;
r2   = p2*mphi*y2/k2;
w    = p2*(1-mphi)*y2/n2;  
dep1 = delm1*km1+delu1*ku1;
dep2 = del2*k2;
dep  = dep1+dep2;
va1  = p1*y1;
va2  = p2*y2;
e2   = (r2-del2-mtauk2)*k2;
ptax = mtauk1*km1;
ctax = mtau1*((rm1-delm1-mtauk1)*km1+ru1*ku1-xu1);
re   = xm1+xu1-delm1*km1;
div  = va1-w*n1-ptax-ctax-xu1-xm1;
e1   = va1-w*n1-delm1*km1-ptax;
atax = mtauc*c+mtaun*(n1+n2)+ctax+ptax+mtau2*e2+mtauk2*k2;
y    = c+xm1+xu1+x2+mg;
ir   = rf;
gnp  = y+mtauc*c-xu1+.04-.611*.04-.0759;
v    = km1f+(1-mtau1)*ku1f;



disp(' ')
disp(' ')
disp(' ')
disp('MEANS OF RATIOS (NOT DIVIDED BY GNP)')
tab1  = [gnp,c,mtauc,xm1,x2,mg,xu1,gnp-y-mtauc*c+xu1, ...
          gnp,dep,gnp-dep,mtauc*c,gnp-dep-mtauc*c, ...
          gnp-dep-mtauc*c,w*(n1+n2),mtauk1*km1+mtauk2*k2,e2,e1, ...
          va1,w*n1,delm1*km1,ptax,e1,ctax,e1-ctax,div,re, ...
          (km1f+ku1f+k2f),km1f,ku1f,k2f, ...
          xm1+xu1+x2,xm1,xu1,x2, ...
          dep,delm1*km1,delu1*ku1,del2*k2, ...
          atax,mtauc*c,mtaun*(n1+n2),ctax,ptax,mtau2*e2,mtauk2*k2, ...
          v]; %/gnp;
tab1  = [tab1,(e1-ctax)*100/v];

str1  = ['  Components of GNP, y+tau_c*c-x_u1-adjustments      %6.3f\n';
         '    Private Consumption less tax, c                  %6.3f\n';
         '    Consumption tax, tau_c*c                         %6.3f\n';
         '    Corporate measured investment, x_m1              %6.3f\n';
         '    Noncorporate investment, x_2                     %6.3f\n';
         '    Government Purchases, g                          %6.3f\n';
         '    LESS: Corporate unmeasured investment, xu1       %6.3f\n';
         '    LESS: Adjustments not made by NIPA               %6.3f\n';
         '                                                          \n';
         '  GNP --> Income                                     %6.3f\n';
         '    - Depreciation, del_m1 *k_m1 + del_2 *k_2        %6.3f\n';
         '    = Net Product                                    %6.3f\n';
         '    - Consumption Taxes, tau_c*c                     %6.3f\n';
         '    = Income                                         %6.3f\n';
         '                                                          \n';
         '  Components of Income, p_1*y_1+p_2*y_2-depreciation %6.3f\n';
         '    Wages, w*(n_1+n_2)                               %6.3f\n';
         '    Property tax, tau_k1 *k_m1+tau_k2 *k_2           %6.3f\n';
         '    Noncorporate Earnings, (r_2-del_2-tau_k2)*k_2    %6.3f\n';
         '    Corporate Earnings (detail below)                %6.3f\n';
         '                                                          \n';
         '  Corporate Output, p_1*y_1                          %6.3f\n';
         '    Corporate wages, w*n_1                           %6.3f\n';
         '    Measured depreciation, del_m1*k_m1               %6.3f\n';
         '    Property tax, tau_k1*k_m1                        %6.3f\n';
         '    Corporate earnings                               %6.3f\n';
         '      Tax, tau_1*[(r_m1-del_m1-tau_k1)*k_m1               \n';
         '                  +r_u1*k_u1-x_u1]                   %6.3f\n';
         '      Earnings after Tax                             %6.3f\n';
         '        Dividends, p_1*y_1-wage-taxes-x_m1-x_u1      %6.3f\n';
         '        Retained earnings, x_m1+x_u1-del_m1*k_m1     %6.3f\n';
         '                                                          \n';

         '  Capital, end of year                               %6.3f\n';
         '    Corporate - measured, k_m1''                      %6.3f\n';
         '    Corporate - unmeasured, k_u1''                    %6.3f\n';
         '    Noncorporate, k_2''                               %6.3f\n';
         '                                                          \n';
         '  Investment                                         %6.3f\n';
         '    Corporate - measured, x_m1                       %6.3f\n';
         '    Corporate - unmeasured, x_u1                     %6.3f\n';
         '    Noncorporate, x_2                                %6.3f\n';
         '                                                          \n';
         '  Depreciation                                       %6.3f\n';
         '    Corporate - measured, del_m1*k_m1                %6.3f\n';
         '    Corporate - unmeasured, del_u1*k_u1              %6.3f\n';
         '    Noncorporate, del_2*k_2                          %6.3f\n';
         '                                                          \n';
         '  Tax Revenues                                       %6.3f\n';
         '    Consumption taxes, tau_c*c                       %6.3f\n';
         '    Wage taxes, tau_n*w*(n_1+n_2)                    %6.3f\n';
         '    Corporate profit taxes, (see formula above)      %6.3f\n';
         '    Corporate property tax, tau_k1*k_m1              %6.3f\n';
         '    Noncorporate taxes, tau_2*e_2                    %6.3f\n';
         '    Noncorporate property tax, tau_k2*k_2            %6.3f\n';
         '                                                          \n';
         '  Value of Corporate Equities                        %6.3f\n';
         '    ==> After-tax earnings/Corporate value *100      %6.3f\n'];

disp(sprintf(str1',tab1))
disp(' ')

[v,rf,w*n1]

