%NX_STD     is the main program for computing equilibria in a stripped-
%           down version of the multi-country economy studied in 
%
%              Technology Capital and the US Current Account
%              by Ellen McGrattan and Ed Prescott
%
%           without technology capital (that is, the standard case).
%           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 Yi = Ai Ki^alpha Hi^(1-alpha), i=u,r
%           Capital accumulates as follows: Ki' = (1-delta) Ki + Xi
%
%           Users input: 
%             * constants for alpha, beta, delta, psi
%             * time series for Ai[t],Ni[t], i=u,r
%             * initial conditions for aggregate Ki[0], i=u,r
%             * 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, 12-14-07

%_____________________________________________________________________________
%
% User inputs
% Benchmark: T=100, alpha=.3, beta=.96, delta=.05, psi=2,
%            o=ones(T,1), Au=o, Ar=.396*o, Nu=o, Nr=8*o,
%            Ku0=1.6049, Kr0=3.4182,
%            iplan=0, B0=0  (or, iplan=1, lam=.78794), imap=1
%_____________________________________________________________________________
%
% 
T     = 100;
o     = ones(T,1);
alpha = .3;
beta  = .96;
delta = .05;
psi   = 2;
Au    = o;
Ar    = .396*ones(T,1); 
Nu    = o;
Nr    = 8*o;          
Ku0   = 1.6049;
Kr0   = 3.4182;
B0    = 0;
lam   = .78974;     
iplan = 0;
imap  = 0;


clc
disp('COMPUTING NET EXPORTS IN MODEL WITHOUT 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('        Yi = Ai Ki^alpha Hi^(1-alpha), i=u,r ')
pause(1)
clc
disp('COMPUTING NET EXPORTS IN MODEL WITHOUT TECHNOLOGY CAPITAL')
disp('---------------------------------------------------------')
disp('    Current parameters: ')
disp(sprintf('       Share parameter, alpha          %8.3f',alpha))
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)))
if iplan==0;
disp(sprintf('       Initial debt of u, Bu[0]        %8.3f',B0))
end;
disp(        '       Initial tangible capital stocks: ')
disp(sprintf('         of US companies, Ku[0]        %8.3f',Ku0))
disp(sprintf('         of ROW companies, Kr[0]       %8.3f',Kr0))
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(' ')
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    = -.0811;
  end;
end;
if alt==2;
  %
  % Add permanent growth in ROW population:
  %
  Nr(43:T)=Nr(42)*ones(T-42,1);
  if adj==1;
    B0    = -.0972;
  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    = .396*Ar;
  if adj==1;
    B0  = 0.1385;
  end;
end;
if alt==4;
  %
  % Add permanent growth in ROW TFP:
  %
  Ar(43:T)=Ar(42)*ones(T-42,1);
  if adj==1;
    B0  = .1669;
  end;
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;
Zu = Au.^(1/(1-alpha));
Zr = Ar.^(1/(1-alpha));
Su = Zu.*Nu;
Sr = Zr.*Nr;
pause(1)

%_____________________________________________________________________________
%
% Call relevant subroutines
%
%_____________________________________________________________________________

h0    = .35;   
Kun   = Ku0*o;
Krn   = Kr0*o;
Hu    = h0*Nu;
Hr    = h0*Nr;
if iplan==1;
  param = [alpha;beta;delta;lam;psi;Ku0;Kr0;T;Au;Ar;Nu;Nr];
  X     = [Kun;Krn;Hu;Hr];
  X     = secant('stdfoc',X,param);
  X1    = zeros(T,4);
  X1(:) = X;
else
  Bn    = B0*o;
  param = [alpha;beta;delta;psi;Ku0;Kr0;B0;T;Au;Ar;Nu;Nr];
  X     = [Kun;Krn;Hu;Hr;Bn];
  X     = secant('decstdfoc',X,param);
  X1    = zeros(T,5);
  X1(:) = X;
  Bn    = X1(:,5);
end;
Kun   = X1(:,1);
Krn   = X1(:,2);
Hu    = X1(:,3);
Hr    = X1(:,4);

%
% Intermediate variables 
%

Nun   = [Nu(2:T);Nu(T)];
Nrn   = [Nr(2:T);Nr(T)];
Ku    = [Ku0;Kun(1:T-1)];
Kr    = [Kr0;Krn(1:T-1)];
Yu    = Au.*Ku.^alpha.*Hu.^(1-alpha);
Yr    = Ar.*Kr.^alpha.*Hr.^(1-alpha);
Cu    = (1-alpha)*Yu.*(Nu-Hu)./(psi*Hu);
Cr    = (1-alpha)*Yr.*(Nr-Hr)./(psi*Hr);
Xu    = Kun-(1-delta)*Ku;
Xr    = Krn-(1-delta)*Kr;
GDPu  = Yu;
GDPr  = Yr;
NXu   = GDPu-Xu-Cu;
NXr   = GDPr-Xr-Cr;
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(NXu./df)-NXu(T)/df(T)*beta/(1-beta);

%
% Plot net exports 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])


%
% Compute B0 = function(lam)
%

if imap==1;
  disp('Computing bond holdings as a function of utility weights ....')
  param = [alpha;beta;delta;lam;psi;Ku0;Kr0;T;Au;Ar;Nu;Nr];
  Lam   = lam-.15:.05:lam+.05;
  lLam  = length(Lam);
  tX    = zeros(T,4);
  XX    = zeros(T*4,lLam);
  i     = find(Lam>=lam1);
  if ~isempty(i);
    i     = i(1);
    X0    = X1(:,1:4);
    for j=i:lLam; 
      param(4) = Lam(j);
      X0       = secant('stdfoc',X0(:),param);
      XX(:,j)  = X0(:); 
    end;
  else
    i   = lLam+1;
  end;
  X0    = X1(:,1:4);
  for j=i-1:-1:1;
    param(4) = Lam(j);
    X0       = secant('stdfoc',X0(:),param);
    XX(:,j)  = X0(:); 
  end;
  for j=1:lLam;
    tX(:) = XX(:,j);
    tKun  = tX(:,1);
    tHu   = tX(:,3);
    tKu   = [Ku0;tKun(1:T-1)];
    tYu   = Au.*tKu.^alpha.*tHu.^(1-alpha);
    tCu   = (1-alpha)*tYu.*(Nu-tHu)./(psi*tHu);
    tXu   = tKun-(1-delta)*tKu;
    tNXu  = tYu-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(tNXu./tdf)-tNXu(T)/tdf(T)*beta/(1-beta);
    By(j) = Bi(j)/tYu(1);
    Yi(j) = tYu(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;    
