%NX_TCAP    is the main program for computing equilibria in a stripped-
%           down version of the multi-country model incorporating 
%           technology capital developed in 
%
%              Technology Capital and the US Current Account
%              by Ellen McGrattan and Ed Prescott
%
%           Here, there are 2 countries (u=US, r=ROW).  
%           Preferences of households in country i are given by
%
%              max sum beta^t {log(Ci[t]/Ni[t])+psi*log(1-Hi[t]/Ni[t])}Ni[t].
%
%           Technologies are Yu = Yud+sigu Yuf, Yr = Yrf+sigr Yrd where
%
%              Yij = Ai (Ni*Mj)^phi (Kij^alpha Hij^(1-alpha))^(1-phi)
%
%           and d=Dell (a US company) and f=Fuji (a non-US company).
%           Capital accumulates as follows: Ki' = (1-delta) Ki + Xi.
%
%           Users input: 
%             * constants for alpha, beta, delta, phi, psi
%             * time series for Ai[t],Ni[t],sigu[t], i=u,r
%             * initial conditions for Kij[0], Mj[0], i=u,r, j=d,f
%             * a choice of solving 
%                 1. a planner's problem with weights u:lam, r:1-lam (iplan=1)
%                 2. a decentralized problem with borrowing and lending and
%                    initial debt levels Bi[0], such that Bu[0]+Br[0]=0
%             * a choice of finding initial B's as functions of utility weights
%                 1. yes, find the functions (imap=1)
%                 2. no, skip it (imap=0)
%
%           A benchmark parameterization is provided.
%

%           Ellen McGrattan, 9-1-07
%           Revised, ERM, 1-17-08

%_____________________________________________________________________________
%
% User inputs
% Benchmark: T=100, alpha=.3, beta=.96, delta=.05, phi=.035, psi=2, o=ones(T,1)
%            Au=o, Ar=.38*o, Nu=o, Nr=8*o, sigu=.75*o, sigr=.75*o,
%            Kud0=1.4303, Kuf0=.00082, Krd0=.00039, Krf0=3.0438,
%            Md0=.1730, Mf0=.3681, iplan=0, B0=0,
%            (or, iplan=1, lam=.789954), imap=0
%_____________________________________________________________________________
%
% 
T     = 100;
o     = ones(T,1);
alpha = .3;
beta  = .96;
delta = .05;
phi   = .035;
psi   = 2.;
Au    = o;
Ar    = .380*o;
Nr    = 8*o;
Nu    = o;
sigu  = .75*o;
sigr  = .75*o;
Kud0  = 1.4303;
Kuf0  = 8.1987e-04;
Krd0  = 3.8526e-04;
Krf0  = 3.0438;
Md0   = 1.7297e-01;
Mf0   = 3.6809e-01;
B0    = 0;
lam   = 0.789954;
iplan = 0;
imap  = 0;
%
% Guesses for per-capita hours---may need to be updated
%
hud0  = 2.9081e-01; 
huf0  = 1.6669e-04;
hrd0  = 3.6824e-05;
hrf0  = 2.9094e-01;

clc
disp('COMPUTING NET EXPORTS IN MODEL WITH TECHNOLOGY CAPITAL')
disp('------------------------------------------------------')
disp('    Countries: u=US, r=ROW')
disp(' ')
disp('    Companies: d=Dell, a US company; f=Fiji, a non-US company')
disp(' ')
disp('    Preferences:  ')
disp('        log(Ci[t]/Ni[t])+psi*log(1-Hi[t]/Ni[t]) ')
disp(' ')
disp('    Technologies: ')
disp('        Yu  = Yud+sigu Yuf ')
disp('        Yr  = Yrf+sigr Yrd ')
disp('        Yij = Ai (Ni Mj)^phi (Kij^alpha Hij^(1-alpha))^(1-phi)')
pause(2)
clc
disp('COMPUTING NET EXPORTS IN MODEL WITH TECHNOLOGY CAPITAL')
disp('------------------------------------------------------')
disp('    Current parameters: ')
disp(sprintf('       Share parameter, alpha          %8.3f',alpha))
disp(sprintf('       Share parameter, phi            %8.3f',phi))
disp(sprintf('       Discount factor, beta           %8.3f',beta))
disp(sprintf('       Depreciation rate, delta        %8.3f',delta))
disp(sprintf('       Utility parameter, psi          %8.3f',psi))
disp(sprintf('       Initial TFP for u, Au[0]        %8.3f',Au(1)))
disp(sprintf('       Initial TFP for r, Ar[0]        %8.3f',Ar(1)))
disp(sprintf('       Initial population of u, Nu[0]  %8.3f',Nu(1)))
disp(sprintf('       Initial population of r, Nr[0]  %8.3f',Nr(1)))
disp(sprintf('       Initial openness of u, sigu[0]  %8.3f',sigu(1)))
disp(sprintf('       Initial openness of r, sigr[0]  %8.3f',sigr(1)))
if iplan==0;
disp(sprintf('       Initial debt of u, Bu[0]        %8.3f',B0))
end;
disp(        '       Initial tangible capital stocks: ')
disp(sprintf('         of Dell in US, Kud[0]         %8.3f',Kud0))
disp(sprintf('         of Fiji in US, Kuf[0]         %8.3f',Kuf0))
disp(sprintf('         of Dell in ROW, Krd[0]        %8.3f',Krd0))
disp(sprintf('         of Fiji in ROW, Krf[0]        %8.3f',Krf0))
disp(        '       Initial technology capital stocks: ')
disp(sprintf('         of Dell, Md[0]                %8.3f',Md0))
disp(sprintf('         of Fiji, Mf[0]                %8.3f',Mf0))
disp(' ')
disp('    Numerical experiments: ')
disp(' ')
disp('       0 = No change in exogenous variables')
disp('       1 = ROW population temporarily rises 18 percent')
disp('       2 = ROW population permanently rises 18 percent')
disp('       3 = ROW TFP temporarily rises 1.2 percent')
disp('       4 = ROW TFP permanently rises 1.2 percent')
disp('       5 = ROW openness parameter increases 20 percent')
disp('       6 = US openness parameter increases 20 percent')
disp('       7 = Both openness parameters increase 20 percent')
disp(' ')
if iplan==0;
disp('    with Bu[0] above or adjusted to avoid initial jumps in investment')
end;
disp(' ')
alt   = input(' Which case? ');
adj   = 0;
if (iplan==0 & alt>0);
  adj = input(' Adjust initial debt? (1=yes) ');
end;

if (alt==1 | alt==2);
  %
  % Add temporary growth in ROW population:
  %
  Nrmax = 1.18;
  Nr    = [linspace(1,Nrmax,38)';
           Nrmax+[.005;.008;.0099;.0103;.0099;.008;.005];
           linspace(Nrmax,1,38)';ones(30,1)];Nr=Nr(1:T);
  Nr    = 8*Nr;
  if adj==1;
    B0    = -.0150;
  end;
end;
if alt==2;
  %
  % Add permanent growth in ROW population:
  %
  Nr(43:T)=Nr(42)*ones(T-42,1);
  if adj==1;
    B0    = -.0173;
  end;
end;
if (alt==3 | alt==4);
  %
  % Add temporary growth in ROW TFP:
  %
  Armax = 1.012;
  Ar    = [linspace(1,Armax,38)';
           Armax+[.00023;.0004;.0005;.00055;.0005;.0004;.00023];
           linspace(Armax,1,38)';ones(30,1)];Ar=Ar(1:T);
  Ar    = .380*Ar;
  if adj==1;
    B0  = .1294;
  end;
end;
if alt==4;
  %
  % Add permanent growth in ROW TFP:
  %
  Ar(43:T)=Ar(42)*ones(T-42,1);
  if adj==1;
    B0  = .1560;
  end;
end;
a     = .75;
b     = .103;
c     = -2;
d     = .05;
if alt==5;
  %
  % Add permanent increase in r's openness parameter
  %
  sigr  = a*(1+b*(tanh(c:d:c+d*(T-1))'-tanh(c)));
  if adj==1;
    B0    = -.0049;
  end;
end;
if alt==6;
  %
  % Add permanent increase in u's openness parameter
  %
  sigu  = a*(1+b*(tanh(c:d:c+d*(T-1))'-tanh(c)));
  if adj==1;
    B0    = .0048;
  end;
end;
if alt==7;
  %
  % Add permanent increase in both openness parameters
  %
  sigu  = a*(1+b*(tanh(c:d:c+d*(T-1))'-tanh(c)));
  sigr  = sigu;
end;

if adj==1;
disp(sprintf('       New initial debt of u, Bu[0]    %8.3f',B0))
end;

t = 1960:1959+T;
if (alt>0 & alt<3);
  figure(1)
  plot(t,Nr./Nu,'-*')
  title('Population of ROW relative to US','Fontsize',14)
end;
if (alt>=3 & alt<=4);
  figure(1)
  plot(t,Ar./Au,'-*')
  title('TFP of ROW relative to US','Fontsize',14)
end;
if (alt>=5)
  figure(1)
  plot(t,[sigr,sigu],'-*')
  title('Openness Parameters','Fontsize',14)
  legend('ROW','US','Location','Northwest')
end;
pause(1)


%_____________________________________________________________________________
%
% Call relevant subroutines
%
%_____________________________________________________________________________
Kudn  = Kud0*o;
Kufn  = Kuf0*o;
Krdn  = Krd0*o;
Krfn  = Krf0*o;
Mdn   = Md0*o;
Mfn   = Mf0*o;
Hud   = hud0*Nu;
Huf   = huf0*Nu;
Hrd   = hrd0*Nr;
Hrf   = hrf0*Nr;
if iplan==1;
  param = [alpha;beta;delta;lam;phi;psi;
           Kud0;Kuf0;Krd0;Krf0;Md0;Mf0;
           T;Au;Ar;Nu;Nr;sigu;sigr];
  X     = [Kudn;Kufn;Krdn;Krfn;Mdn;Mfn;Hud;Huf;Hrd;Hrf];
  X     = secant('tcapfoc',X,param);
  X1    = zeros(T,10);
  X1(:) = X;
else
  Bn    = B0*o;
  param = [alpha;beta;delta;phi;psi;
           Kud0;Kuf0;Krd0;Krf0;Md0;Mf0;B0;
           T;Au;Ar;Nu;Nr;sigu;sigr];
  X     = [Kudn;Kufn;Krdn;Krfn;Mdn;Mfn;Hud;Huf;Hrd;Hrf;Bn];
  X     = secant('dectcapfoc',X,param);
  X1    = zeros(T,11);
  X1(:) = X;
  Bn    = X1(:,11);
end;
Kudn  = X1(:,1);
Kufn  = X1(:,2);
Krdn  = X1(:,3);
Krfn  = X1(:,4);
Mdn   = X1(:,5);
Mfn   = X1(:,6);
Hud   = X1(:,7);
Huf   = X1(:,8);
Hrd   = X1(:,9);
Hrf   = X1(:,10);

%
% Intermediate variables 
%

Nun   = [Nu(2:T);Nu(T)];
Nrn   = [Nr(2:T);Nr(T)];
Kud   = [Kud0;Kudn(1:T-1)];
Kuf   = [Kuf0;Kufn(1:T-1)];
Krd   = [Krd0;Krdn(1:T-1)];
Krf   = [Krf0;Krfn(1:T-1)];
Md    = [Md0;Mdn(1:T-1)];
Mf    = [Mf0;Mfn(1:T-1)];
Kun   = Kudn+Kufn;
Ku    = Kud+Kuf;
Krn   = Krdn+Krfn;
Kr    = Krd+Krf;
Xud   = Kudn-(1-delta)*Kud;
Xuf   = Kufn-(1-delta)*Kuf;
Xrd   = Krdn-(1-delta)*Krd;
Xrf   = Krfn-(1-delta)*Krf;
Xd    = Mdn-(1-delta)*Md;
Xf    = Mfn-(1-delta)*Mf;
Xu    = Xud+Xuf;
Xr    = Xrd+Xrf;
Hu    = Hud+Huf;
Hr    = Hrd+Hrf;
Yud   = Au.*(Nu.*Md).^phi.*(Kud.^alpha.*Hud.^(1-alpha)).^(1-phi);
Yuf   = Au.*sigu.*(Nu.*Mf).^phi.*(Kuf.^alpha.*Huf.^(1-alpha)).^(1-phi);
Yrd   = Ar.*sigr.*(Nr.*Md).^phi.*(Krd.^alpha.*Hrd.^(1-alpha)).^(1-phi);
Yrf   = Ar.*(Nr.*Mf).^phi.*(Krf.^alpha.*Hrf.^(1-alpha)).^(1-phi);
Cu    = (1-alpha)*(1-phi)*Yud.*(Nu-Hu)./(psi*Hud);
Cr    = (1-alpha)*(1-phi)*Yrf.*(Nr-Hr)./(psi*Hrf);
GDPu  = Yud+Yuf-Xd;
GDPr  = Yrd+Yrf-Xf;
NFIu  = phi*(Yrd-Yuf);
NFIr  = phi*(Yuf-Yrd);
GNPu  = GDPu+NFIu;
GNPr  = GDPr+NFIr;
NXu   = GDPu-Cu-Xu;
NXr   = GDPr-Cr-Xr;
CAu   = GNPu-Cu-Xu;
CAr   = GNPr-Cr-Xr;
lam1  = Cu(1)/Nu(1)/(Cr(1)/Nr(1)+Cu(1)/Nu(1));
Cun   = [Cu(2:T);Cu(T)];
cu    = Cu./Nu;
cun   = Cun./Nun;
rbn   = (cun./cu/beta)-1;
rb    = [1/beta-1;rbn(1:T-1)];
df    = cumprod(1+rb);
B1    = -sum(CAu./df)-CAu(T)/df(T)*beta/(1-beta);

%
% Plot net exports and current account relative to GDP
% (assuming period 1 is 1960)
%

figure(1)
load data.dat
plot(t,NXu./GDPu*100,'b-','Linewidth',2)
line(data(:,1),data(:,2),'Color','b','Marker','o')
title('Net Exports Relative to GDP (%)','Fontsize',14)
legend('Theory','US Data')
ylabel('Percent')
axis([1960,2010,-8,4])

figure(2)
plot(t,GDPu./Nu*Nu(1)/GDPu(1)*100,'Color','b','Linewidth',2)
line(t,GDPr./Nr*Nr(1)/GDPr(1)*100,'Color','r','Linewidth',2)
title('Model Predictions for US and ROW Per Capita GDP (1960=100)', ...
      'Fontsize',14)
legend('US','ROW')
ylabel('Index')
axis([1960,2010,90,105])

figure(3)
prdu = GDPu./Hu;
prdr = GDPr./Hr;
plot(t,prdu/prdu(1)*100,'Color','b','Linewidth',2)
line(t,prdr/prdr(1)*100,'Color','r','Linewidth',2)
title('Model Predictions for US and ROW Labor Productivity (1960=100)', ...
      'Fontsize',14)
legend('US','ROW')
ylabel('Index')
axis([1960,2010,99,104])

figure(4)
cuy = Cu./GDPu*100;
cry = Cr./GDPr*100;
plot(t,cuy,'Color','b','Linewidth',2)
line(t,cry,'Color','r','Linewidth',2)
title('Model Predictions for US and ROW Consumption Shares (%)', ...   
      'Fontsize',14)
legend('US','ROW')
ylabel('Percent')
axis([1960,2010,75,95])

figure(5)
xuy = Xu./GDPu*100;
xry = Xr./GDPr*100;
plot(t,xuy,'Color','b','Linewidth',2)
line(t,xry,'Color','r','Linewidth',2)
title('Model Predictions for US and ROW Tangible Investment Shares (%)', ...   
      'Fontsize',14)
legend('US','ROW')
ylabel('Percent')
axis([1960,2010,5,25])

figure(6)
mdy = Md./GDPu;
mfy = Mf./GDPr;
plot(t,mdy/mdy(1)*100,'Color','b','Linewidth',2)
line(t,mfy/mfy(1)*100,'Color','r','Linewidth',2)
title(['Model Predictions for Technology Capital Stocks Relative ', ...
       'to own country GDP (1960=100)'],'Fontsize',14)
legend('US Companies','ROW Companies')
ylabel('Index')
axis([1960,2010,90,110])

%
% Compute B0 = function(lam)
%

if imap==1;
  disp('Computing bond holdings as a function of utility weights ....')
  param = [alpha;beta;delta;lam;phi;psi;Kud0;Kuf0;Krd0;Krf0;Md0;Mf0;T;
           Au;Ar;Nu;Nr;sigu;sigr];
  Lam   = lam-.15:.05:lam+.05;
  lLam  = length(Lam);
  tX    = zeros(T,10);
  XX    = zeros(T*10,lLam);
  i     = find(Lam>=lam1);
  if ~isempty(i); 
    i   = i(1);
    X0    = X1(:,1:10);
    for j=i:lLam; 
      param(4) = Lam(j);
      X0       = secant('tcapfoc',X0(:),param);
      XX(:,j)  = X0(:); 
    end;
  else
    i   = lLam+1;
  end;
  X0    = X1(:,1:10);
  for j=i-1:-1:1;
    param(4) = Lam(j);
    X0       = secant('tcapfoc',X0(:),param);
    XX(:,j)  = X0(:); 
  end;
  for j=1:lLam;
    tX(:) = XX(:,j);
    tKudn = tX(:,1);
    tKufn = tX(:,2);
    tKrdn = tX(:,3);
    tMdn  = tX(:,5);
    tMfn  = tX(:,6);
    tHud  = tX(:,7);
    tHuf  = tX(:,8);
    tHrd  = tX(:,9);
    tKud  = [Kud0;tKudn(1:T-1)];
    tKuf  = [Kuf0;tKufn(1:T-1)];
    tKrd  = [Krd0;tKrdn(1:T-1)];
    tMd   = [Md0;tMdn(1:T-1)];
    tMf   = [Mf0;tMfn(1:T-1)];
    tKun  = tKudn+tKufn;
    tKu   = tKud+tKuf;
    tXu   = tKun-(1-delta)*tKu;
    tXd   = tMdn-(1-delta)*tMd;
    tHu   = tHud+tHuf;
    tYud  = Au.*(Nu.*tMd).^phi.*(tKud.^alpha.*tHud.^(1-alpha)).^(1-phi);
    tYuf  = Au.*sigu.*(Nu.*tMf).^phi.*(tKuf.^alpha.*tHuf.^(1-alpha)).^(1-phi);
    tYrd  = Ar.*sigr.*(Nr.*tMd).^phi.*(tKrd.^alpha.*tHrd.^(1-alpha)).^(1-phi);
    tCu   = (1-alpha)*(1-phi)*tYud.*(Nu-tHu)./(psi*tHud);
    tGDPu = tYud+tYuf-tXd;
    tGNPu = tGDPu+phi*(tYrd-tYuf);
    tCAu  = tGNPu-tCu-tXu;
    tXu   = tKun-(1-delta)*tKu;
    tNXu  = tGDPu-tXu-tCu;
    tCun  = [tCu(2:T);tCu(T)];
    tcu   = tCu./Nu;
    tcun  = tCun./Nun;
    trbn  = (tcun./tcu/beta)-1;
    trb   = [1/beta-1;trbn(1:T-1)];
    tdf   = cumprod(1+trb);
    Bi(j) = -sum(tCAu./tdf)-tCAu(T)/tdf(T)*beta/(1-beta);
    By(j) = Bi(j)/tGDPu(1);
    Yi(j) = tGDPu(1);
  end;

%
% Plot lam versus B0
%

figure(2)
plot(Lam,Bi,'-*')
title('Relation of Initial Bond Holdings and Utility Weight','Fontsize',14)
xlabel('Weight on US Utility in Planner''s Problem','Fontsize',14)
ylabel('Initial Bond Position of US','Fontsize',14)

figure(3)
plot(Lam,By*100,'-*')
title('Relation of Initial Debt Levels and Utility Weight','Fontsize',14)
xlabel('Weight on US Utility in Planner''s Problem','Fontsize',14)
ylabel('Initial Bond Position of US Relative to US GDP (%)','Fontsize',14)
end;    
