%BENCHMARK  computes the transition path for the benchmark model in:
%           
%              "Why are Married Women Working So Much?"
%               by Larry Jones, Rodolfo Manuelli, and Ellen McGrattan
%
%           as well as special cases which can be run by uncommenting 
%           relevant lines at the bottom of the input file (benchmarki).
%

%
% Iterate from benchmark cases:
% 

bparam=[5.0000;20.0000; 0.3000; 0.2000; 0.9039; 0.3409; 0.3409; 0.3409; 
        0.1905; 0.1905; 0.1858; 0.1858; 1.0000; 1.0000; 1.0000; 1.0000; 
        0.3265; 0.2200; 1.0000; 0.2200; 1.0000; 2.0000; 0.7673;-0.2867; 
        0.5859; 4.0000; 1.0000; 1.0000; 1.0000; 0.2000; 0.5000; 0.4500; 
        0.2300; 0.4500; 0.2300; 0.6000; 0.2000; 0.6000; 0.2000; 0.0000;
        0.0000; 1.0000; 1.0000; 1.0000; 1.0000; 1.0000; 1.0000];

bvars= [0.0973; 0.4914; 0.3451; 0.3760; 0.2500; 0.0819; 0.0811; 0.0175; 
        0.0293; 0.0175; 0.0294; 0.0106; 0.0178; 0.0422; 0.0211; 0.0211; 
        0.0081; 0.0030; 0.0054; 0.0074; 0.0230; 0.0106; 0.0204; 0.0074; 
        0.0230; 0.0106; 0.0204; 0.0011; 0.0022; 0.5000];

%
% ... to the desired case:
% 
benchmarki
bparam(22) = util;
bparam(46) = res1;
bparam(47) = res2;
param      = [nyr;nper;np;nfs;beta;delh;dele;delk;kap1;zet1;kap2;zet2;
              b1f;b2f;b1m;b2m;lam;theta1;mtech;theta2;htech;util;psi1;
              psi2;psi3;sigma;gpar1;gpar2;gpar3;g;tauk;taulfp;taulmp;
              taulfs;taulms;Tp;Tfs;Tdp;Tdfs;gama;gamp;qk1;qk2;qh;qe;
              res1;res2];
vars       = [l1fp; l1mp; l1fs; l1ms; l2fp; l2fs; l2ms;
              c1fp; c1mp; c1fs; c1ms; c2fp; c2mp;
              k1p;  k1fs; k1ms; k2p;  k2fs; k2ms;
              hfp;  hmp;  hfs;  hms;  efp;  emp;  efs;  ems; T; Td; relw];

npar0      = 47;
nx0        = 30;
nx1        = 29;
npar1      = 45;
npcomp     = length(pcomp);
xon        = ones(nx0,1);
pon        = zeros(npar0,1);
if ~isempty(pcomp);
  xon(vnocomp) = zeros(npcomp,1);
  pon(pcomp)   = ones(npcomp,1);
end;
SS         = [bvars(find(xon));bparam(find(pon))];

BPARAMS    = [xon;bvars;pon;bparam];
PARAMS     = [xon;vars;pon;param];
nwgt = 10;
nwgt = 500;
wgt  = linspace(1,0,nwgt);
for i=1:nwgt;
  SS   = secant('ssfoc',SS,BPARAMS*wgt(i)+PARAMS*(1-wgt(i)));
end;
if any(imag(SS));
  SS   = secant('ssfoc',real(SS),PARAMS);
end;

disp('NEW PARAMETER VALUES:')
STR=[ ...
 '    nyr    = %g;\n'
 '    nper   = %g;\n'
 '    np     = %g;\n'
 '    nfs    = %g;\n'         
 '    beta   = %g;\n'   
 '    delh1  = %g;\n'
 '    delh2  = %g;\n'
 '    delk   = %g;\n' 
 '    kap1   = %g;\n'
 '    zet1   = %g;\n'
 '    kap2   = %g;\n'
 '    zet2   = %g;\n'
 '    b1f    = %g;\n'   
 '    b2f    = %g;\n'        
 '    b1m    = %g;\n'       
 '    b2m    = %g;\n'      
 '    lam    = %g;\n'
 '    theta1 = %g;\n'
 '    mtech  = %g;\n'
 '    theta2 = %g;\n'
 '    htech  = %g;\n'
 '    util   = %g;\n'
 '    psi1   = %g;\n'
 '    psi2   = %g;\n'
 '    psi3   = %g;\n'
 '    sigma  = %g;\n'
 '    gpar1  = %g;\n'
 '    gpar2  = %g;\n'   
 '    gpar3  = %g;\n'  
 '    g      = %g;\n'
 '    tauk   = %g;\n'
 '    taulfp = %g;\n'
 '    taulmp = %g;\n'
 '    taulfs = %g;\n'
 '    taulms = %g;\n'
 '    Tp     = %g;\n'
 '    Tfs    = %g;\n'
 '    Tdp    = %g;\n'
 '    Tdfs   = %g;\n'
 '    gama   = %g;\n'
 '    gamp   = %g;\n'
 '    qk1    = %g;\n'
 '    qk2    = %g;\n'
 '    qh     = %g;\n'
 '    qe     = %g;\n'];

vars(find(xon))  = SS(1:length(find(xon)));
param(find(pon)) = SS(length(find(xon))+1:nx0);
if param(46) == 1;
  param(10) = param(9);
end;
if param(47) == 1;
  param(12) = param(11);
end;

disp(sprintf(STR',param(1:npar1)))
%((1+param(40))^(1-(1-param(23))*(1-param(26)))/param(5))^(1/nyr)
if ~any(imag(SS));
  bvars=vars;
  bparam=param;
else
  disp('WARNING: imaginary numbers found in steady state calculation.')
end;
%save junk bvars bparam
%return
VARS   = [  '    Married Female Market Hours           %7.4f\n';
            '    Married Male Market Hours             %7.4f\n';
            '    Single Female Market Hours            %7.4f\n';
            '    Single Male Market Hours              %7.4f\n';
            '    Married Female Home Hours             %7.4f\n';
            '    Single Female Home Hours              %7.4f\n';
            '    Single Male Home Hours                %7.4f\n';
            '    Married Female Market Consumption     %7.4f\n';
            '    Married Male Market Consumption       %7.4f\n';
            '    Single Female Market Consumption      %7.4f\n';
            '    Single Male Market Consumption        %7.4f\n';
            '    Married Female Home Consumption       %7.4f\n';
            '    Married Male Home Consumption         %7.4f\n';
            '    Married Capital in the Market         %7.4f\n';
            '    Single Female Market Capital          %7.4f\n';
            '    Single Male Market Capital            %7.4f\n';
            '    Married Homes Capital                 %7.4f\n';
            '    Single Female Home Capital            %7.4f\n';
            '    Single Male Home Capital              %7.4f\n';
            '    Married Female Human Capital Type h   %7.4f\n';
            '    Married Male Human Capital Type h     %7.4f\n';
            '    Single Female Human Capital Type h    %7.4f\n';
            '    Single Male Human Capital Type h      %7.4f\n';
            '    Married Female Human Capital Type eta %7.4f\n';
            '    Married Male Human Capital Type eta   %7.4f\n';
            '    Single Female Human Capital Type eta  %7.4f\n';
            '    Single Male Human Capital Type eta    %7.4f\n';
            '    Government Transfers                  %7.4f\n';
            '    Discrimination Transfers              %7.4f\n';
            '    Relative wage                         %7.4f\n'];
disp(sprintf(VARS',vars))
%return

Param    = [param(1:npar1);Mtech;Htech;Taulfp;Taulmp;Taulfs;Taulms;
            Qk1;Qk2;Qh;Qe;csfp;csmp;csfs;csms;tsfp;tsmp;tsfs;tsms];

%
% Initial guess for transition path: sequence of steady states.
%

nper    = Param(2);
X       = zeros(nx1,nper);
X(:,1)  = vars(1:nx1);
for i=2:nper;
  Param(19)  = Param(npar1+0*nper+i);
  Param(21)  = Param(npar1+1*nper+i);
  Param(32)  = Param(npar1+2*nper+i);
  Param(33)  = Param(npar1+3*nper+i);
  Param(34)  = Param(npar1+4*nper+i);
  Param(35)  = Param(npar1+5*nper+i);
  Param(42)  = Param(npar1+6*nper+i);
  Param(43)  = Param(npar1+7*nper+i);
  Param(44)  = Param(npar1+8*nper+i);
  Param(45)  = Param(npar1+9*nper+i);
  X(:,i)     = secant('ssfoc2',X(:,i-1),Param);
end;
disp('Steady state calculations complete. ')

X      = [X;zeros(4,nper)];
nx1    = nx1+4;

%
% Compute transition path.
%

disp('Computing the transition path ....  ')
cap0   = X(14:27,1);
Y      = ssecant('trfoc',X(:),[Param;cap0]);
if abs(sum(imag(Y)))>0;
  Y    = ssecant('trfoc',real(Y),[Param;cap0]);
end;
tem    = zeros(nx1,nper);
tem(:) = Y;
Y      = tem;

domore=0;
if domore==1;
nmore  = 8;
Y      = [Y,Y(:,nper)*ones(1,nmore)];
Mtech  = [Mtech;Mtech(nper)*ones(nmore,1)];
Htech  = [Htech;Htech(nper)*ones(nmore,1)];
Taulfp = [Taulfp;Taulfp(nper)*ones(nmore,1)];
Taulmp = [Taulmp;Taulmp(nper)*ones(nmore,1)];
Taulfs = [Taulfs;Taulfs(nper)*ones(nmore,1)];
Taulms = [Taulms;Taulms(nper)*ones(nmore,1)];
Qk1    = [Qk1;Qk1(nper)*ones(nmore,1)];
Qk2    = [Qk2;Qk2(nper)*ones(nmore,1)];
Qh     = [Qh;Qh(nper)*ones(nmore,1)];
Qe     = [Qe;Qe(nper)*ones(nmore,1)];
nper   = nper+nmore;
Param  = [param(1:npar1);Mtech;Htech;Taulfp;Taulmp;Taulfs;Taulms;
          Qk1;Qk2;Qh;Qe;csfp;csmp;csfs;csms;tsfp;tsmp;tsfs;tsms];
Param(2) = nper;
Y      = ssecant('trfoc',Y(:),[Param;cap0]);
tem    = zeros(nx1,nper);
tem(:) = Y;
Y      = tem;
end;

%
% Report results.
%



%
% Solution:
%

l1fp     = Y(1,:)';
l1mp     = Y(2,:)';
l1fs     = Y(3,:)';
l1ms     = Y(4,:)';
l2fp     = Y(5,:)';
l2fs     = Y(6,:)';
l2ms     = Y(7,:)';
c1fp     = Y(8,:)';
c1mp     = Y(9,:)';
c1fs     = Y(10,:)';
c1ms     = Y(11,:)';
c2fp     = Y(12,:)';
c2mp     = Y(13,:)';
k1pn     = Y(14,:)';
k1fsn    = Y(15,:)';
k1msn    = Y(16,:)';
k2pn     = Y(17,:)';
k2fsn    = Y(18,:)';
k2msn    = Y(19,:)';
hfpn     = Y(20,:)';
hmpn     = Y(21,:)';
hfsn     = Y(22,:)';
hmsn     = Y(23,:)';
efpn     = Y(24,:)';
empn     = Y(25,:)';
efsn     = Y(26,:)';
emsn     = Y(27,:)';
T        = Y(28,:)';
Td       = Y(29,:)';
mufp     = Y(30,:)';
mump     = Y(31,:)';
mufs     = Y(32,:)';
mums     = Y(33,:)';

%
% Parameters used to compute other variables:
%

nyr      = Param(1);
nper     = Param(2);
np       = Param(3);
nfs      = Param(4);
nms      = 1-2*np-nfs;
beta     = Param(5);
delh1    = Param(6);
delh2    = Param(7);
delk     = Param(8);
kap1     = Param(9);
zet1     = Param(10);
kap2     = Param(11);
zet2     = Param(12);
b1f      = Param(13);
b2f      = Param(14);
b1m      = Param(15);
b2mf     = Param(16);
lam      = Param(17);
theta1   = Param(18);
theta2   = Param(20);
util     = Param(22);
psi1     = Param(23);
psi2     = Param(24);
psi3     = Param(25);
sigma    = Param(26);
gpar1    = Param(27);
gpar2    = Param(28);
gpar3    = Param(29);
g        = Param(30);
tauk     = Param(31);
Tp       = Param(36)*T;
Tfs      = Param(37)*T;
Tms      = T-Tp-Tfs;
Tdp      = Param(38)*Td;
Tdfs     = Param(39)*Td;
Tdms     = Td-Tdp-Tdfs;
gama     = Param(40);
gamp     = Param(41);
Mtech    = Param(46+0*nper:45+1*nper);
Htech    = Param(46+1*nper:45+2*nper);
Taulfp   = Param(46+2*nper:45+3*nper);
Taulmp   = Param(46+3*nper:45+4*nper);
Taulfs   = Param(46+4*nper:45+5*nper);
Taulms   = Param(46+5*nper:45+6*nper);
Qk1      = Param(46+6*nper:45+7*nper);
Qk2      = Param(46+7*nper:45+8*nper);
Qh       = Param(46+8*nper:45+9*nper);
Qe       = Param(46+9*nper:45+10*nper);
grate    = ((1+gama)*(1+gamp)).^[0:nper-1]';
csfp     = csfp(1)*ones(nper,1)./grate;
csmp     = csmp(1)*ones(nper,1)./grate;
csfs     = csfs(1)*ones(nper,1)./grate;
csms     = csms(1)*ones(nper,1)./grate;

%
% Other variables:
%

k1p      = [cap0(1); k1pn(1:nper-1)];
k1fs     = [cap0(2); k1fsn(1:nper-1)];
k1ms     = [cap0(3); k1msn(1:nper-1)];
k2p      = [cap0(4); k2pn(1:nper-1)];
k2fs     = [cap0(5); k2fsn(1:nper-1)];
k2ms     = [cap0(6); k2msn(1:nper-1)];
hfp      = [cap0(7); hfpn(1:nper-1)];
hmp      = [cap0(8); hmpn(1:nper-1)];
hfs      = [cap0(9); hfsn(1:nper-1)];
hms      = [cap0(10);hmsn(1:nper-1)];
efp      = [cap0(11);efpn(1:nper-1)];
emp      = [cap0(12);empn(1:nper-1)];
efs      = [cap0(13);efsn(1:nper-1)];
ems      = [cap0(14);emsn(1:nper-1)];
z1fp     = b1f*hfp.^kap1.*efp.^zet1.*l1fp;   
z1mp     = b1m*hmp.^kap1.*emp.^zet1.*l1mp;   
z1fs     = b1f*hfs.^kap1.*efs.^zet1.*l1fs;   
z1ms     = b1m*hms.^kap1.*ems.^zet1.*l1ms;   
z2fp     = b2f*hfp.^kap2.*efp.^zet2.*l2fp;   
z2fs     = b2f*hfs.^kap2.*efs.^zet2.*l2fs;   
z2ms     = b2m*hms.^kap2.*ems.^zet2.*l2ms;   
kbar     = nfs*k1fs+nms*k1ms+np*k1p;
zbar     = nfs*z1fs+nms*z1ms+np*z1fp+np*z1mp;
[F1,w,r]            = prodfn(zbar,kbar,theta1,Mtech);
[c1p,F2zp,F2kp]     = prodfn(z2fp,k2p, theta2,Htech);
[c2fs,F2zfs,F2kfs]  = prodfn(z2fs,k2fs,theta2,Htech);
[c2ms,F2zms,F2kms]  = prodfn(z2ms,k2ms,theta2,Htech);
xk1p     = (1+gama)*(1+gamp)*k1pn -(1-delk)*k1p;
xk1fs    = (1+gama)*(1+gamp)*k1fsn-(1-delk)*k1fs;
xk1ms    = (1+gama)*(1+gamp)*k1msn-(1-delk)*k1ms;
xk2p     = (1+gama)*(1+gamp)*k2pn -(1-delk)*k2p;
xk2fs    = (1+gama)*(1+gamp)*k2fsn-(1-delk)*k2fs;
xk2ms    = (1+gama)*(1+gamp)*k2msn-(1-delk)*k2ms;
xhfp     = (1+gama)*(1+gamp)*hfpn-(1-delh)*hfp;
xefp     = (1+gama)*(1+gamp)*efpn-(1-dele)*efp;
xhmp     = (1+gama)*(1+gamp)*hmpn-(1-delh)*hmp;
xemp     = (1+gama)*(1+gamp)*empn-(1-dele)*emp;
xhfs     = (1+gama)*(1+gamp)*hfsn-(1-delh)*hfs;
xefs     = (1+gama)*(1+gamp)*efsn-(1-dele)*efs;
xhms     = (1+gama)*(1+gamp)*hmsn-(1-delh)*hms;
xems     = (1+gama)*(1+gamp)*emsn-(1-dele)*ems;
qshare   = Qk2.*(nfs*xk2fs+nms*xk2ms+np*xk2p)./F1;
l1avg    = nfs*l1fs+nms*l1ms+np*(l1mp+l1fp);


Tax      = np *[tauk*r.*k1p-tauk*delk*k1p+Taulfp.*w.*z1fp+Taulmp.*w.*z1mp]+ ...
           nfs*[tauk*r.*k1fs-tauk*delk*k1fs+Taulfs.*w.*z1fs]+ ...
           nms*[tauk*r.*k1ms-tauk*delk*k1ms+Taulms.*w.*z1ms];
Dtax     = np *(Taulfp-Taulms).*w.*z1fp + nfs *(Taulfs-Taulms).*w.*z1fs- ...
           np *(Taulms-Taulmp).*w.*z1mp;
Gtax     = Tax-Dtax;
Gov      = g*F1;
Tran     = Tax-Gov;
Td       = Dtax;
T        = Gtax-Gov;

%
% DATA:
%


% Hours for married females & males, single females and males, ages 25-64
hours = [ ...
    1950    8.0024   41.6833   28.8010   31.4252
    1960   10.8532   40.1376   28.0608   29.1782
    1970   13.7832   39.5603   27.0204   28.0432
    1980   18.0191   38.6713   28.2560   29.9582
    1990   23.3182   39.5396   29.5931   31.1863
    2000   24.7334   38.7731   29.2903   30.9101
    2005   24.7388   38.6083   29.1703   32.2942];

% Gender gap constructed by Blau (1998)
blau =[ 1969, .562;
        1979, .583;
        1989, .682;
        1994, .717;
        2009, .736];

% Home investment share, 1950=100
xk2divy=[ ...
    1955  100.00
    1965   93.73
    1975   98.83
    1985   92.63
    1995   90.29
    2005   92.96];

%
% FIGURES:
%
DH   = 84;   % Discretionary hours = 84 hours per week
Y0   = 1930;
Y0   = 1950;
figure(1)
plot(Y0:nyr:Y0-1+nyr*nper,l1fp*DH,'b-',hours(:,1),hours(:,2),'b+', ...
     Y0:nyr:Y0-1+nyr*nper,l1mp*DH,'r-',hours(:,1),hours(:,3),'r+', ...
     Y0:nyr:Y0-1+nyr*nper,l1ms*DH,'g-',hours(:,1),hours(:,5),'g+', ...
     Y0:nyr:Y0-1+nyr*nper,l1fs*DH,'y-',hours(:,1),hours(:,4),'y+')
title(['Market Hours of Work ', ...
       '(Model = Solid, Data = +, MM = r, MF = b, SF = y, SM = g)'])
text(2005,40,'Married Males')
text(2005,30,'Singles')
text(2005,18,'Married Females')

figure(2)
plot(Y0:nyr:Y0-1+nyr*nper,(1-Taulfp).*z1fp.*l1mp./ ...
                           ((1-Taulmp).*z1mp.*l1fp),'b-', ...
     Y0:nyr:Y0-1+nyr*nper,(1-Taulfs).*z1fs.*l1ms./ ...
                           ((1-Taulms).*z1ms.*l1fs),'r-', ...
     Y0:nyr:Y0-1+nyr*nper,((1-Taulfp).*z1fp./l1fp.*np+ ...
                             (1-Taulfs).*z1fs./l1fs.*nfs)./ ...
                            ((1-Taulmp).*z1mp./l1mp.*np+ ...
                             (1-Taulms).*z1ms./l1ms.*nms),'g-', ...
     blau(:,1),blau(:,2),'k*')
Relw = ((1-Taulfp).*z1fp./l1fp*np+ ...
        (1-Taulfs).*z1fs./l1fs*nfs)./ ...
       ((1-Taulmp).*z1mp./l1mp*np+ ...
       (1-Taulms).*z1ms./l1ms*nms);
Relwf = ((1-Taulfs).*z1fs./l1fs)./ ...
        ((1-Taulfp).*z1fp./l1fp);
Relwm = ((1-Taulms).*z1ms./l1ms)./ ...
        ((1-Taulmp).*z1mp./l1mp);

title(['Relative wages rates (Married=blue, Single=red, All=green)'])

figure(3)
plot(Y0:nyr:Y0-1+nyr*nper,z1fp.*l1mp./(z1mp.*l1fp),'b-', ...
     Y0:nyr:Y0-1+nyr*nper,z1fs.*l1ms./(z1ms.*l1fs),'r-', ...
     Y0:nyr:Y0-1+nyr*nper,(z1fp./l1fp.*np+z1fs./l1fs.*nfs)./ ...
                            (z1mp./l1mp.*np+z1ms./l1ms.*nms),'g-', ...
     blau(:,1),blau(:,2),'g+')
title(['Relative pre-tax wages rates (Married  =blue, Single=red, All=green)'])

figure(4)
plot(Y0:nyr:Y0-1+nyr*nper,(l1fp+l2fp)*DH,'b-', ...
     Y0:nyr:Y0-1+nyr*nper, l1mp*DH,'r-', ...
     Y0:nyr:Y0-1+nyr*nper,(l1ms+l2ms)*DH,'g-', ...
     Y0:nyr:Y0-1+nyr*nper,(l1fs+l2fs)*DH,'y-')
title('Total Hours of Work (MM = r, MF = b, SF = y, SM = g)')

figure(5)
orient tall
subplot(311)
plot(Y0:nyr:Y0-1+nyr*nper,l1fp*DH,'b-',hours(:,1),hours(:,2),'b+', ...
     Y0:nyr:Y0-1+nyr*nper,l1mp*DH,'r-',hours(:,1),hours(:,3),'r+', ...
     Y0:nyr:Y0-1+nyr*nper,l1ms*DH,'g-',hours(:,1),hours(:,5),'g+', ...
     Y0:nyr:Y0-1+nyr*nper,l1fs*DH,'y-',hours(:,1),hours(:,4),'y+')
title(['Market Hours of Work ', ...
       '(Model = Solid, Data = +, MM = r, MF = b, SF = y, SM = g)'])
axis([1950,2010,5,45])
subplot(312)
plot(Y0:nyr:Y0-1+nyr*nper,((1-Taulfp).*z1fp./l1fp.*np+ ...
                             (1-Taulfs).*z1fs./l1fs.*nfs)./ ...
                            ((1-Taulmp).*z1mp./l1mp.*np+ ...
                             (1-Taulms).*z1ms./l1ms.*nms),'b-', ...
                             blau(:,1),blau(:,2),'b+')
title(['Relative wage rate'])
axis([1950,2010,.4,1])
subplot(313)
plot(Y0:nyr:Y0-1+nyr*nper,qshare/qshare(1)*100,'b-', ...
     xk2divy(:,1),xk2divy(:,2), 'b+')
title(['Share of Home Investment Expenditures (1950=100)'])
axis([1950,2010,50,150])

figure(6)
orient tall
subplot(311)
plot(Y0:nyr:Y0-1+nyr*nper,hfp/hmp(1),'b-', ...
     Y0:nyr:Y0-1+nyr*nper,hmp/hmp(1),'r-', ...
     Y0:nyr:Y0-1+nyr*nper,hms/hmp(1),'g-', ...
     Y0:nyr:Y0-1+nyr*nper,hfs/hmp(1),'y-');
title('Human Capital Relative to MM in 1950')
axis([1950,2010,0,1])
subplot(312)
plot(Y0:nyr:Y0-1+nyr*nper,Relwf,'b', ...
     Y0:nyr:Y0-1+nyr*nper,Relwm,'r')
title('Relative single to married wage rates (M = r, F = b)')
axis([1950,2010,.9,1.2])
subplot(313)
plot(Y0:nyr:Y0-1+nyr*nper,F1./l1avg *l1avg(1)/F1(1)*100)
title('Average productivity relative to 2% trend (F1/avg l) with 1950=100')
axis([1950,2010,90,110])

disp(' ')
disp(' ')
disp('Line for Table3.tex:')
disp('--------------------')
tem= [[l1fp(1),l1mp(1),l1fs(1),l1ms(1)]*DH,Relw(1)*100,     ...
 [l1fp(11),l1mp(11),l1fs(11),l1ms(11)]*DH,Relw(11)*100,qshare(11)/qshare(1)*100];
disp(sprintf(' &&%5.1f&&%5.1f&&%5.1f&&%5.1f&&%5.0f&&%5.1f&&%5.1f&&%5.1f&&%5.1f&&%5.0f&&%5.0f\\cr',tem'))


disp(' ')
disp(' ')
disp('Original Line for Table3.tex (for 1990):')
disp('----------------------------------------')
tem= [[l1fp(1),l1mp(1),l1fs(1),l1ms(1)]*DH,Relw(1)*100,     ...
 [l1fp(9),l1mp(9),l1fs(9),l1ms(9)]*DH,Relw(9)*100,qshare(9)/qshare(1)*100];
disp(sprintf(' &&%5.1f&&%5.1f&&%5.1f&&%5.1f&&%5.0f&&%5.1f&&%5.1f&&%5.1f&&%5.1f&&%5.0f&&%5.0f\\cr',tem'))
