%
% FIXEDPAR   loads data from IOtables and stores fixed parameters in 
%            file fixedpar.dat
%


%
% Shares from IO table
% (Run setupio.m first)
%
load ../IOtable/omega.dat
load ../IOtable/gamma.dat
load ../IOtable/zeta.dat
load ../IOtable/nu.dat
load ../IOtable/capshr.dat
load ../IOtable/cftratio.dat


%
% Distribute weight on "Gc" goods proportionally to other commodities
%

omega(15)  = .001;
omega      = omega/sum(omega);

%
% Adjust capital shares to mimic ratios of (pivec.*prodxi)./(ptvec.*prodxt):
% 
%    pivec(s)*prodxi(s) = sum_l q(l)*XI(l,s) = sum of col s of CFTI
%    ptvec(s)*prodxt(s) = sum_l p(l)*XT(l,s) = sum of col s of CFTT
%

phi        =[0.0200;0.0420;0.1000;0.3280;0.5470; ...
             0.5410;0.3740;0.1400;0.5440;0.0800; ...
             0.7330;0.2950;0.3200;0.2870;0.2600].*capshr';
theta      = capshr'-phi;

%
% Remaining parameters
%

n          = 15;
alpha      = 1.0001;
beta       = 0.995;
cbar       = 1;
deli       = 0.008*1*ones(n,1);
delt       = 0.008*1*ones(n,1);
gn         = 0.0025*1;
gz         = 0.005*1;
p1         = omega(1);
p1         = 1;
psi        = 1.20;
rho        = 1.0001;
z1         = 1.00*ones(n,1);
z2         = 1.00*ones(n,1);
g          = 1e-15*ones(n,1);
tauc       = 0.065;
taud       = 0.144;
tauh       = 0.382;
tauk       = 0.003*1;
taup       = 0.33;
taux       = 0.00;
z          = 1.00;

xguess = [ ...
   0.001291367911315
   0.108655124515771
   0.399352988436393
   0.273353644664319
   0.510713126128129
   0.178564586863085
   0.190914845646467
   0.547743553045316
   0.219893767207328
   0.295395314929304
   0.099755452868877
   0.177440193519153
   0.308465113090407
   0.203815908125713
   0.134958606683901
   0.000147321096024
   0.001831675469171
   0.003769568014064
   0.001946230682825
   0.136896475028933
   0.029132275882732
   0.005619380354396
   0.003031131425181
   0.348339563169064
   0.004254806428413
   0.755061852602025
   0.009390847769797
   0.000534397646483
   0.012599843182992
   0.044962317299706
   0.001444773719911
   0.002327120894790
   0.002854239319852
   0.012271057324424
   0.033599057204334
   0.017167065105336
   0.016663986928303
   0.009280124054806
   0.007962513908102
   0.025127527685714
   0.047741201749724
   0.030603022421931
   0.011751479729589
   0.009887357004470
   0.038934069881604
];

param         = [n;alpha;beta;cbar;deli;delt;g;gamma(:);gn;gz;nu(:);omega;
                 p1;phi;psi;rho;tauc;taud;tauh;tauk;taup;taux;theta;z1;z2;
                 zeta(:);z];

x1            = ss_fullmodel(xguess,param);

%
% Steady state:
%

o             = ones(n,1);
betg          = beta*(1+gz)^(-alpha);
gam           = sum(gamma)';
lshare        = 1-theta-phi-gam;
thet1         = 1-theta;
aggcnrm       = x1(1);
pvec          = [p1;x1(2:n)];
h2h1          = x1(n+1:2*n);
h             = x1(2*n+1:3*n);
pmvec         = gam.*prod(((pvec*o')./(gamma+1e-15)).^(gamma./(o*gam')))';
ptvec         = prod(((pvec*o')./(zeta+1e-15)).^zeta)';
aggp          = sum(omega.^rho.*pvec.^(1-rho))^(1/(1-rho));
cnrm          = (aggp*omega./pvec).^rho * aggcnrm;
c             = cnrm*cbar;
aggc          = (sum(omega.*c.^((rho-1)/rho)))^(rho/(rho-1));
aggl          = 1-sum(h);
m1y           = gam.*pvec./pmvec;
h1            = h./(1+h2h1);
h2            = h-h1;
qxipy         = h2h1;
xiy           = qxipy.^(1-phi).*(z2./z1).^lshare;
qp            = qxipy./xiy;
qvec          = qp.*pvec;
pivec         = prod(((qvec*o')./(nu+1e-15)).^nu)';
kiy           = betg*phi.*(qvec.*xiy+pvec)./((1-betg*(1-deli)).*pivec);
ret           = ((1+taux)*ptvec./pvec-betg*(1-(1-taup)*(delt+tauk) + ...
                (1-delt).*((1+taux)*ptvec./pvec-1)))/(betg*(1-taup));
k1ty          = theta./ret;
y             = (k1ty.^theta.*kiy.^phi.*m1y.^gam).^(o./lshare).*z1.*h1*z;
k1t           = k1ty.*y;
k2t           = h2h1.*k1t;
kt            = k1t+k2t;
m1            = m1y.*y;
xi            = xiy.*y;
ki            = kiy.*y;
qxi           = qvec.*xi;
qxip          = qxi./pvec;
prodxt        = ((1+gz)*(1+gn)-1+delt).*kt;
prodxi        = ((1+gz)*(1+gn)-1+deli).*ki;
m1mat         = gamma.*((o./pvec)*pvec').*(o*y');
m2mat         = gamma.*((o./pvec)*pvec').*(o*qxip');
xtmat         = zeta.*((o./pvec)*ptvec').*(o*prodxt');
ximat         = nu.*((o./qvec)*pivec').*(o*prodxi');
Sbar          = [log(z1);log(z2);log(g);tauc;taud;tauh;tauk;taup;taux;log(z)];
cftmodel      = (pivec.*prodxi)./(ptvec.*prodxt);
fixedpar      = [param;xguess;pvec;qvec;y;xi;1-aggl]; 

isave         = input('Overwrite results in fixedpar.dat? (1=yes,0=no) ');
if (isave==1)
  save fixedpar.dat fixedpar -ascii
end;

return

%
% When iterating on parameters, use:
%

clc
disp('Parameters are stored in vector fixedpar')
disp('Adjust phi/theta until ratios below are aligned')
disp(' ')
disp(' ')
disp(' Intangible to Tangible Ratio')
disp('=============================')
disp(' Industry   Data      Model  ')
disp('-----------------------------')
for i=1:n;
  disp(sprintf('%5g    %7.2f   %7.2f',i,cftratio(i),cftmodel(i))); 
end;
disp('=============================')
disp('NOTE: pivec(s)*prodxi(s)       ')
disp('         = sum_l q(l)*XI(l,s)  ')
disp('         = sum of col s of CFTI')
disp('      ptvec(s)*prodxt(s)       ')
disp('         = sum_l p(l)*XT(l,s)  ')
disp('         = sum of col s of CFTE+CFTS ')

