%ED     this program computes certain steady state results found by running
%       GMOD_UC.m.  See the latter for details of the model.
%         

%       Ellen McGrattan, 7-13-00

iout  = input('Output file name (e.g., ed17.out): ','s');
fid   = fopen(iout,'w');

iexp  = 0;
istop = 0;
while istop==0;

iexp  = iexp+1;
load ed.in;
Dgamma   = ed(1)/4;           
Deta     = ed(2)/4;          
Dc       = ed(3);           
Dxm1     = ed(4);           
Dx2      = ed(5);           
Dg       = ed(6);           
Dw1      = ed(7);           
Dw2      = ed(8);           
Dr1x     = ed(9);           
Dr2      = ed(10);           
Dkm1     = ed(11)*4;        
Dk2      = ed(12)*4;       
Dibt     = ed(13);          
Dprof    = ed(14);          
Dtxprof  = ed(15);          
n        = ed(16);           
rho      = ed(17);
sigma    = ed(18);
mtau2    = ed(19);
mtaun    = ed(20);

% Back out model parameters:

gamma    = Dgamma;
eta      = Deta;
mtauc    = Dibt/Dc;
mtau1    = Dtxprof/Dprof;
delm1    = Dxm1/Dkm1 - (gamma+eta+gamma*eta);
delu1    = delm1;
del2     = Dx2/Dk2 - (gamma+eta+gamma*eta);
mphi     = Dr2/(Dw2+Dr2);
Mr2      = Dr2/Dk2;
beta     = (1+gamma)^sigma/((1-mtau2)*(Mr2-del2)+1);
betg     = beta*(1+gamma)^(-sigma);
Mru1     = (1-mtau2)*(Mr2-del2)+delu1;
Mrm1     = (Mru1-delu1)/(1-mtau1)+delm1;
alpha    = (Dw1+Dw2+Mrm1*Dkm1+Dr2-Dc-Dxm1-Dx2-Dg)/(Dxm1-Mru1*Dkm1);
Dxu1     = Dxm1*alpha;
Dy       = Dc+Dxm1+Dxu1+Dx2+Dg;
Dr1      = Dr1x+Dxu1;
Dva1     = Dr1+Dw1;
Dva2     = Dr2+Dw2;
Dku1     = Dkm1*alpha;
mthetm   = Mrm1*Dkm1/Dva1;
mthetu   = Mru1*Dku1/Dva1;
Mw       = ((1-mphi)*Dva2+(1-mthetm-mthetu)*Dva1)/n;
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;


% Steady state values of interest:

My       = Mp1*My1+Mp2*My2;      
Mgnp     = (My+mtauc*Dc-Dxu1)/1.046;
Mv       = (1+grate)*(1+alpha)*Dkm1 /4;
Mi       = (1+gamma)^sigma/beta-1;
Mprof    = (Mrm1-delm1)*Dkm1+(Mru1-delu1)*Dku1;


% Print out results:
disp(' ')
disp(' ')
disp(sprintf('EXPERIMENT %g ',iexp))
disp(' ')
disp(' ')
disp('Inputs Used (Annualized): ')
disp(sprintf('  gamma                =  %5.3f',4*Dgamma)) 
disp(sprintf('  eta                  =  %5.3f',4*Deta)) 
disp(sprintf('  c                    =  %5.3f',Dc)) 
disp(sprintf('  xm1                  =  %5.3f',Dxm1)) 
disp(sprintf('  x2                   =  %5.3f',Dx2)) 
disp(sprintf('  g                    =  %5.3f',Dg)) 
disp(sprintf('  w*n1                 =  %5.3f',Dw1)) 
disp(sprintf('  w*n2                 =  %5.3f',Dw2)) 
disp(sprintf('  rm1*km1+ru1*ku1-xu1  =  %5.3f',Dr1x)) 
disp(sprintf('  r2*k2                =  %5.3f',Dr2)) 
disp(sprintf('  km1                  =  %5.3f',Dkm1/4))
disp(sprintf('  k2                   =  %5.3f',Dk2/4))
disp(sprintf('  ibt                  =  %5.3f',Dibt))
disp(sprintf('  corp. prof.          =  %5.3f',Dprof))
disp(sprintf('  tax on corp. prof.   =  %5.3f',Dtxprof))
disp(sprintf('  time at work         =  %5.3f',n))
disp(sprintf('  rho                  =%7.3f',rho))
disp(sprintf('  sigma                =  %5.3f',sigma))
disp(sprintf('  tau_2                =  %5.3f',mtau2))
disp(sprintf('  tau_n                =  %5.3f',mtaun))
disp(' ')
disp(' ')
disp('Model Parameters Backed Out: ')
disp(sprintf('  mtauc=ibt/c          =  %5.3f',mtauc))
disp(sprintf('  mtau1=txprof/prof    =  %5.3f',mtau1))
disp(sprintf('  delm1 from x,k,growth=  %5.3f',delm1*4))
disp(sprintf('  delu1=delm1          =  %5.3f',delu1*4))
disp(sprintf('  del2 from x,k,growth =  %5.3f',del2*4))
disp(sprintf('  theta2 from R2,VA2   =  %5.3f',mphi))
disp(sprintf('  r2 from R2,K2        =  %5.3f',Mr2))
disp(sprintf('  beta from i=return_2 =  %5.3f',beta))
disp(sprintf('  ru1 from i=return_u1 =  %5.3f',Mru1))
disp(sprintf('  rm1 from i=return_m1 =  %5.3f',Mrm1))
disp(sprintf('  alpha from prod=inc  =  %5.3f',alpha))
disp(sprintf('  xu1=alpha*xm1        =  %5.3f',Dxu1))
disp(sprintf('  theta_m1 from RM1,VA1=  %5.3f',mthetm))
disp(sprintf('  theta_u1 from RU1,VA1=  %5.3f',mthetu))
disp(sprintf('  w from (W1+W2)/n     =  %5.3f',Mw))
disp(sprintf('  n1 from W1,W2,n      =  %5.3f',Mn1))
disp(sprintf('  n2 from W1,W2,n      =  %5.3f',Mn2))
disp(sprintf('  psi from UL/UC=w-tax =  %5.3f',psi))
disp(sprintf('  a1 from firm FOC     =  %5.3f',a1))
disp(sprintf('  a2 from firm FOC     =  %5.3f',a2))
disp(' ')
disp(' ')
disp('Steady States of Interest: ')
disp(sprintf('  interest rate        =  %5.3f',Mi*400))
disp(sprintf('  gnp                  =  %5.3f',Mgnp))
disp(sprintf('  v/gnp                =  %5.3f',Mv/Mgnp))
disp(sprintf('  corporate profits/gnp=  %5.3f',Mprof/Mgnp))
disp(sprintf('  depreciation_m1/gnp  =  %5.3f',delm1*Dkm1/Mgnp))
disp(sprintf('  depreciation_u1/gnp  =  %5.3f',delu1*Dku1/Mgnp))
disp(sprintf('  depreciation_2/gnp   =  %5.3f',del2*Dk2/Mgnp))

sv  = input('Save results? (y or n) ','s');

if (sv=='y' | sv=='Y');

  fprintf(fid,' \n');
  fprintf(fid,' \n');
  fprintf(fid,'EXPERIMENT %g \n',iexp);
  fprintf(fid,' \n');
  fprintf(fid,' \n');
  fprintf(fid,'Inputs Used (Annualized): \n');
  fprintf(fid,'  gamma                =  %5.3f\n',4*Dgamma);
  fprintf(fid,'  eta                  =  %5.3f\n',4*Deta);
  fprintf(fid,'  c                    =  %5.3f\n',Dc);
  fprintf(fid,'  xm1                  =  %5.3f\n',Dxm1); 
  fprintf(fid,'  x2                   =  %5.3f\n',Dx2); 
  fprintf(fid,'  g                    =  %5.3f\n',Dg); 
  fprintf(fid,'  w*n1                 =  %5.3f\n',Dw1); 
  fprintf(fid,'  w*n2                 =  %5.3f\n',Dw2); 
  fprintf(fid,'  rm1*km1+ru1*ku1-xu1  =  %5.3f\n',Dr1x); 
  fprintf(fid,'  r2*k2                =  %5.3f\n',Dr2); 
  fprintf(fid,'  km1                  =  %5.3f\n',Dkm1/4);
  fprintf(fid,'  k2                   =  %5.3f\n',Dk2/4);
  fprintf(fid,'  ibt                  =  %5.3f\n',Dibt);
  fprintf(fid,'  corp. prof.          =  %5.3f\n',Dprof);
  fprintf(fid,'  tax on corp. prof.   =  %5.3f\n',Dtxprof);
  fprintf(fid,'  time at work         =  %5.3f\n',n);
  fprintf(fid,'  rho                  =%7.3f\n',rho);
  fprintf(fid,'  sigma                =  %5.3f\n',sigma);
  fprintf(fid,'  tau_2                =  %5.3f\n',mtau2);
  fprintf(fid,'  tau_n                =  %5.3f\n',mtaun);
  fprintf(fid,' \n');
  fprintf(fid,' \n');
  fprintf(fid,'Model Parameters Backed Out: \n');
  fprintf(fid,'  mtauc=ibt/c          =  %5.3f\n',mtauc);
  fprintf(fid,'  mtau1=txprof/prof    =  %5.3f\n',mtau1);
  fprintf(fid,'  delm1 from x,k,growth=  %5.3f\n',delm1*4);
  fprintf(fid,'  delu1=delm1          =  %5.3f\n',delu1*4);
  fprintf(fid,'  del2 from x,k,growth =  %5.3f\n',del2*4);
  fprintf(fid,'  theta2 from R2,VA2   =  %5.3f\n',mphi);
  fprintf(fid,'  r2 from R2,K2        =  %5.3f\n',Mr2);
  fprintf(fid,'  beta from i=return_2 =  %5.3f\n',beta);
  fprintf(fid,'  ru1 from i=return_u1 =  %5.3f\n',Mru1);
  fprintf(fid,'  rm1 from i=return_m1 =  %5.3f\n',Mrm1);
  fprintf(fid,'  alpha from prod=inc  =  %5.3f\n',alpha);
  fprintf(fid,'  xu1=alpha*xm1        =  %5.3f\n',Dxu1);
  fprintf(fid,'  theta_m1 from RM1,VA1=  %5.3f\n',mthetm);
  fprintf(fid,'  theta_u1 from RU1,VA1=  %5.3f\n',mthetu);
  fprintf(fid,'  w from (W1+W2)/n     =  %5.3f\n',Mw);
  fprintf(fid,'  n1 from W1,W2,n      =  %5.3f\n',Mn1);
  fprintf(fid,'  n2 from W1,W2,n      =  %5.3f\n',Mn2);
  fprintf(fid,'  psi from UL/UC=w-tax =  %5.3f\n',psi);
  fprintf(fid,'  a1 from firm FOC     =  %5.3f\n',a1);
  fprintf(fid,'  a2 from firm FOC     =  %5.3f\n',a2);
  fprintf(fid,' \n');
  fprintf(fid,' \n');
  fprintf(fid,'Steady States of Interest: \n');
  fprintf(fid,'  interest rate        =  %5.3f\n',Mi*400);
  fprintf(fid,'  gnp                  =  %5.3f\n',Mgnp);
  fprintf(fid,'  v/gnp                =  %5.3f\n',Mv/Mgnp);
  fprintf(fid,'  corporate profits/gnp=  %5.3f\n',Mprof/Mgnp);
  fprintf(fid,'  depreciation_m1/gnp  =  %5.3f\n',delm1*Dkm1/Mgnp);
  fprintf(fid,'  depreciation_u1/gnp  =  %5.3f\n',delu1*Dku1/Mgnp);
  fprintf(fid,'  depreciation_2/gnp   =  %5.3f\n',del2*Dk2/Mgnp);

end;

cn  = input('Continue experimenting? (y or n) ','s');

if (cn=='n' | cn=='N');
  istop = 1;
end;

end;

fclose('all');







