
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file solves the continuous time model of Fanelli Straub (2021) %
%     " A Theory of Foreign Exchange Interventions "                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all

% header
width = '3in';
height = '2in';
scale  = 0.8; % slides: 0.65. Paper: 0.8
suffix = '_paper'; % change to '' for slides!
globalstyle = 'very thick';

% paths
H.root_path = '';
H.path_tikz = [H.root_path 'tikz/'];
H.compile_pdf_folder     = 'tikz_pdfs/';
H.compile_pdf_tikzfolder = 'tikz/';
H.pdflatexpath = '/usr/local/texlive/2018/bin/x86_64-darwin'; % location of pdflatex
fprintf(['\n\n *** NOTE *** \n Location of pdflatex set to: ' H.pdflatexpath '\n If this location is incorrect the tikz and pdf production will fail.\n Please set write_tikz and compile_pdf to 0 or adjust pdflatexpath.\n'])

% Set these two options to 1 to produce tikz files and pdfs
%  NOTE: this *only* works if above path 'pdflatexpath' was set to the
%  location of pdflatex! Otherwise, choose these to be 0
compile_pdf = 1;
write_tikz  = 1;


%% Set up
P.rho      = 0.075;   % discount rate
P.delta    = 0;       % rate of decay in bonds = frequency of adjustment (not used)
P.yNT      = 1;       % total NT endowment (only important for the price)
P.alpha    = 0.051148; % openness parameter in preferences
P.mu       = 0.2;     % population share of poor
P.yT       = 0.234;   % total T endowment
P.yT_extra = 1.7764e-16; % not used
P.c_lbar   = 0.19271; % subsistence level of tradables
chi        = 0.078;   % share of income going to poor
P.gamma_T  = chi;     % share of T income earned by poor
P.gamma_NT = chi;     % share of NT income earned by poor
P.sigma    = 2;       % inverse EIS

% set model solution options
P.mode      = 'opt'; % could be 'opt','eqm','opt_iniER', 'tau0', 'opt_res_acc', 'opt_unc', 'opt_LC'
P.break_indeterminacy = 'res_const'; % how to break asymmetry if tau = 0 ? 'res_const' sets reserves constant during tau = 0 stretches, 'res_linear' interpolates reserves linearly
P.addl_vars = 'StoneGeary'; % pick 'StoneGeary' or 'homebias' for what model to use to compute additional outcome variables
P.res_lbar  = -Inf; % lower bound for reserves (unused)
P.res_lbar_rstardot = @(t) 0; % important for lower bound: what is rstar dot as fn of t


%%% Objective and constraint
% auxiliary objects
P.yT_h2m  = P.gamma_T*P.yT/P.mu;    % per capita T income of poor 
P.yNT_h2m = P.gamma_NT*P.yNT/P.mu;  % per capita NT income of poor

% set parameters to compute objective
a0 = 1-P.mu; b0 = P.mu; c0 = 1-P.alpha;
b1 = P.mu*P.alpha/(P.gamma_NT*(1-P.mu)*(1-P.alpha))*(P.c_lbar-P.yT_h2m);
c1 = b1*(1-P.alpha)*P.gamma_NT;
h0 = -P.mu*P.alpha*(P.c_lbar-P.yT_h2m)/(1-(1-P.alpha)*P.gamma_NT) - P.yT - P.yT_extra + P.c_lbar;
h1 = 1-P.mu + P.gamma_NT*(1-P.mu)*(1-P.alpha)/(1-(1-P.alpha)*P.gamma_NT);

% objective V
syms a
V_sym = a0*a + b0*log(exp(a)-b1) - c0*log(exp(a)-c1);
if P.sigma ~= 1
    V_sym = (exp(a)-c1)^(c0*(P.sigma-1)) * ( ...
                a0*exp(a*(1-P.sigma)) + b0 * (a0*c1/b0/b1/(1-c1/b1))^(1-P.sigma) * (exp(a)-b1)^(1-P.sigma) ...
                )/(1-P.sigma) - 1/(1-P.sigma);
end
P.V   = matlabFunction(V_sym);
P.Vp  = matlabFunction(diff(V_sym));
P.Vpp = matlabFunction(diff(V_sym,2));

% cost of consumption function h
syms a
h_sym = h0 + h1*exp(a);
P.h   = matlabFunction(h_sym);
P.hp  = matlabFunction(diff(h_sym));
P.hpp = matlabFunction(diff(h_sym,2));
P.h_inv = matlabFunction(finverse(h_sym));

% RER function + its inverse
%  - important if one wishes to set initial exchange rate targets
pp0 = (1-P.alpha)/(P.yNT*(1-(1-P.alpha)*P.gamma_NT));
pp1 = (1-P.mu)/P.alpha;
pp2 = P.mu*(P.c_lbar-P.yT_h2m);
P.p_fn     = @(x) pp0 * ( pp1*x - pp2); % price of NT
P.pp_fn    = @(x) pp0 * pp1;            % derivative w.r.t. total spending x 
P.p_inv_fn = @(p) (p/pp0 + pp2) / pp1;
P.RER_fn   = @(a) 1./( (1/P.alpha)^P.alpha * (P.p_fn(exp(a))/(1-P.alpha)).^(1-P.alpha) );
P.RER_inv_fn = @(rer) log( P.p_inv_fn(  (1/rer*P.alpha^P.alpha*(1-P.alpha).^(1-P.alpha))^(1/(1-P.alpha))  ) );


if b1 < 0
    warning 'Need c_lbar > yT_h2m to have the right sign of intervention.'
    P.yT_h2m
end

% initialize
previousrun = 0; k = 1;

% common parameters:
P.abstol = 1e-9; P.reltol = 1e-9;
T = 25;
nfa0 = 0;
Gamma = 9;
Gamma_list = [1 9 50];
mass = 0.004;        % of home intermediaries
GammaH_ratio = 1/3;  % Ratio GammaH / Gamma
shock = struct('rstar', @(t) P.rho - 0.02*(t<=5).*(1-t/5));  % set common rstar shock
P.D   = create_demand(Gamma,Gamma*GammaH_ratio,0,mass);
P_common = P;

Z_common = struct;
Z_common.fields_to_plot = {'nfaGDP','costGDP','RERss','bI','resGDP',...
                            'tau','rstar','bIH','privnfaGDP', 'z_R_ss', 'z_h2m_ss','p_ss'}; %'zeta'
Z_common.field_labels   = {'Net foreign assets (NFA) / GDP','Carry trade cost / GDP','Real exchange rate','bI','Reserves / GDP',...
                            'Interest rate spread $\\tau_t$','World interest rate $r^*$','bIH','NFA excl. reserves / GDP', 'Ricardian expenditure $z_t$', 'Poor expenditure $z^P_t$', 'Nontradable price'};
Z_common.stepsize       = [0,0,0,0,0,0,1,0,0,0,0,0,0,0]; % if zero, then Matlab's stepsize is used
Z_common.yaxis_percent  = ones(1,12);
Z_common.x_label        = 'years';
Z_common.first_is_norm  = [0,0,0,0,0,0,0,0,0,0,0,0]; % use first tdlist as normalization? if = 2, first entry is not shown
Z_common.max_T          = min(10,T); % highest T shown
Z_common.height         = height;
Z_common.width          = width;
Z_common.globalstyle    = globalstyle;
Z_common.compile_pdf    = compile_pdf;
Z_common.write_tikz     = write_tikz;
Z_common.scale          = scale;

Z_debug = Z_common;
Z_debug.fields_to_plot = {'nfaGDP','costGDP','RERss','bI','resGDP','tau','rstar','bIH','privnfaGDP','zeta','mu'};
Z_debug.field_labels   = {'Net foreign assets (NFA) / GDP','Carry trade cost / GDP','Real exchange rate','bI','Reserves / GDP','Interest rate spread $\\tau_t$','World interest rate $r^*$','bIH','NFA excl. reserves / GDP','zeta','mu'};
Z_debug.stepsize       = [0,0,0,0,0,0,1,0,0,0,0,0,0];
Z_debug.first_is_norm  = [0,0,0,0,0,0,0,0,0,0,0];
Z_debug.yaxis_percent  = ones(1,11);


%% Figure 1: optimal intervention time series
clearvars tdlist
P = P_common;
k=1; P.mode = 'eqm'; % first run equilibrium
P.D   = create_demand(Gamma,Gamma*GammaH_ratio,0,mass);
tdlist{k} = compute_td(nfa0,T,P,shock,0);k=k+1;

P.mode = 'opt'; % optimal for various gammas
for Gamma_here = Gamma_list
    P.D   = create_demand(Gamma_here,Gamma_here*GammaH_ratio,0,mass);
    tdlist{k} = compute_td(nfa0,T,P,shock,tdlist{k-1});k=k+1;
end

%%% plot it
Z = Z_common;
Z.plot_name      = 'opt';
Z.colorlist      = {'black!70','green!20','green!50','green'};
Z.stylelist      = {'solid','solid','solid','solid'};
Z.legend_        = {'No intervention',['$\\Gamma = ' num2str(Gamma_list(1)) '$'],['$\\Gamma = ' num2str(Gamma_list(2)) '$'],['$\\Gamma = ' num2str(Gamma_list(3)) '$']};
Z.plot_name = [Z.plot_name suffix]; plot_this(Z,tdlist,H);


%% Figure 3: PV of costs for smooth ER and optimum
P = P_common; P.abstol = 1e-7; P.reltol = 1e-7; P.iniER = 0;
Gamma_lb = 0.01; Gamma_ub = 100; fn = @(x) log(x); fn_inv = @(x) exp(x);
Gamma_list_here = fn_inv(linspace(fn(Gamma_ub),fn(Gamma_lb),30));
cost.opt = zeros(1,length(Gamma_list_here));
cost.opt_iniER = zeros(1,length(Gamma_list_here));
quasifiscalcost.opt = zeros(1,length(Gamma_list_here));
quasifiscalcost.opt_iniER = zeros(1,length(Gamma_list_here));
modes = fields(cost);

% compute
for i = 1:2
    k = 1; P.mode = modes{i}; clearvars tdlist
    for Gamma_here = Gamma_list_here
        P.D   = create_demand(Gamma_here,Gamma_here*GammaH_ratio,0,mass);
        if k > 1
            tdlist{k} = compute_td(nfa0,T,P,shock,tdlist{k-1});k=k+1;
        else
            tdlist{k} = compute_td(nfa0,T,P,shock,0);k=k+1;
        end
        % compute costs
        discount = (1+P.rho).^(-tdlist{k-1}.t);
        PV_cost  = trapz( tdlist{k-1}.t, tdlist{k-1}.costGDP .* discount );
        PV_quasifiscalcost = trapz( tdlist{k-1}.t, tdlist{k-1}.resGDP .* tdlist{k-1}.tau .* discount );
        cost.(modes{i})(k-1) = PV_cost * 100;
        quasifiscalcost.(modes{i})(k-1) = PV_quasifiscalcost * 100;
    end
end

O        = struct;
O.scale  = scale;
O.data   = {[Gamma_list_here', cost.(modes{1})'], [Gamma_list_here', cost.(modes{2})'], ...
            [Gamma_list_here', quasifiscalcost.(modes{1})'], [Gamma_list_here', quasifiscalcost.(modes{2})']};
O.plot_name = ['cost'];
O.xaxis  = {Gamma_lb, Gamma_ub};
O.yaxis  = {0.00011, max(quasifiscalcost.(modes{2}))};
O.title  = 'Intervention cost (present value, \\% of GDP)';
O.colors = {'green','red','green','red'};
O.styles = {'solid','solid','densely dashed','densely dashed'};
O.legend = {'Optimal','Smooth exchange rate','Optimal (quasi-fiscal cost)', 'Smooth exchange rate (quasi-fiscal cost)'};
O.legend_cols = 2; % max 2 legend columns
O.xlabel = 'financial friction $\\Gamma$';
O.globalstyle = 'line width = 1pt';
O.width  = '4in';
O.height = '3in';
% O.width_off  = 1;
% O.height_off = 1;
O.yaxis_percent = 1;
%     O.xaxis_percent = '';
O.xtick = {0.0001,0.001,0.01,0.1,1,10,100};
O.ytick = {0.0001,0.001,0.01,0.1,1,10,100}; %create_ticks(O.yaxis{1},O.yaxis{2},stepsize_cost,true);
O.compile_pdf    = compile_pdf; 
O.compile_pdf_folder = H.compile_pdf_folder;
O.compile_pdf_paperpath = H.root_path;
O.compile_pdf_tikzfolder = H.compile_pdf_tikzfolder;
O.axis = 'loglogaxis';

write_tikz_file(O, H.path_tikz);


%% Figure 2: smooth ER vs unconstrained optimal
P = P_common; clearvars tdlist
k=1; P.mode = 'opt';
tdlist{k} = compute_td(nfa0,T,P,shock,previousrun);k=k+1;
P.mode = 'opt_iniER'; P.iniER = 0;
tdlist{k} = compute_td(nfa0,T,P,shock,tdlist{k-1});k=k+1;

%%% plot this
Z = Z_common;
Z.plot_name      = 'smooth_ER';
Z.colorlist      = {'green','red'};
Z.stylelist      = {'solid','solid'};
Z.legend_        = {['Optimal ($\\Gamma = ' num2str(Gamma) '$)'],'Smooth exchange rate'};
Z.plot_name = [Z.plot_name suffix]; plot_this(Z,tdlist,H);


%% Figure 4: forward guidance puzzle
P = P_common; P.abstol = 1e-7; P.reltol = 1e-7;
P.D = create_demand(Gamma,Gamma*GammaH_ratio,0,0); % exclude mass from this edxercise, since it's differenced out in regressions!
P.res_goal   = 5*0.008859; % goal to accumulate reserves after time T (if 'opt_res_acc'). < 1% here because we want to have 1% of *GDP*
P.res_goal_T = 1;
P.mode = 'eqm'; k = 1; clearvars tdlist
for i = [0,3,6]
    shock_here = struct('reserves',@(t) (t>=i) * (t<=i+P.res_goal_T) * P.res_goal );
    tdlist{k} = compute_td(nfa0,T,P,shock_here,previousrun);k=k+1;
end

%%% plot this
Z = Z_common;
Z.plot_name      = 'fwd_guid_puzzle';
Z.colorlist      = {'green','green!50','green!20'};
Z.stylelist      = {'solid','solid','solid'};
Z.legend_        = {'Intervention at t = 0','Intervention at t = 3','Intervention at t = 6'};
Z.max_T          = 10;
Z.height         = '2.5in';
Z.width          = '3.5in';
Z.plot_name = [Z.plot_name suffix]; plot_this(Z,tdlist,H);


%% Figure 5: Limited commitment
clearvars tdlist
P = P_common; legend = {}; P.D   = create_demand(Gamma,Gamma*GammaH_ratio,0,mass);
k=1; P.mode = 'eqm'; % as reference
tdlist{k} = compute_td(nfa0,T,P,shock,0);k=k+1;
legend{k-1} = 'No commitment ($=$ no intervention)';

P.mode = 'opt'; % run full commitment optimum
tdlist{k} = compute_td(nfa0,T,P,shock,0);k=k+1;
legend{k-1} = 'Full commitment';

P.mode = 'opt_LC'; % limited commitment
LC_horizons = [2.5 1 0.1];
for LC_horizon = LC_horizons
    P.LC_dates = [LC_horizon:LC_horizon:5];
    tdlist{k} = compute_td(nfa0,T,P,shock,0);k=k+1;
    legend{k-1} = ['Horizon = ' num2str(LC_horizon)];
end

P.mode = 'eqm';
tdlist{k} = tdlist{1}; k=k+1;
legend{k-1} = 'Horizon = 0';

%%% plot it
Z = Z_common;
Z.first_is_norm(3) = 1; % normalize RER relative to eqm
Z.field_labels{3} = 'Effect of intervention on the real exchange rate';
Z.plot_name      = 'opt_LC';
Z.colorlist      = {'black','green','blue!50!white!50!green','blue!50!white','blue!50!white!50!black','black'};
Z.stylelist      = {'solid','solid','solid','solid','solid','solid'};
Z.legend_        = legend;
Z.height         = '2.5in';
Z.width          = '3.5in';
Z.plot_name = [Z.plot_name suffix]; plot_this(Z,tdlist,H);


%% Figure 8: Simple FXI experiment swap to calibrate Gamma
P = P_common; P.abstol = 1e-7; P.reltol = 1e-7; Gamma = 9;
P.D = create_demand(Gamma,Gamma*GammaH_ratio,0,0); % exclude mass from this edxercise, since it's differenced out in regressions!
P.res_goal   = 1.0524/100; % goal to accumulate reserves after time T (if 'opt_res_acc'). < 1% here because we want to have 1% of *GDP*
P.res_goal_T = 1;220/365;  % what time the reserve goal must be achieved by
shock_here = struct('reserves',@(t) (t<=P.res_goal_T) * P.res_goal );
P.mode = 'eqm'; k = 1; clearvars tdlist
tdlist{k} = compute_td(nfa0,T,P,shock_here,previousrun);k=k+1;

%%% plot this
Z = Z_common;
Z.plot_name      = 'res_acc';
Z.colorlist      = {'green'};
Z.stylelist      = {'solid'};
Z.legend_        = {'Simple reserve swap'};
Z.max_T          = 2;
Z.plot_name = [Z.plot_name suffix]; plot_this(Z,tdlist,H);


%% Figure 10: Learning over time (declining Gamma)
clearvars tdlist
P = P_common;
k=1; P.mode = 'eqm'; % first run equilibrium
P.D   = create_demand(Gamma,Gamma*GammaH_ratio,0,mass);
tdlist{k} = compute_td(nfa0,T,P,shock,0);k=k+1;

P.mode = 'opt'; % optimal for various gammas
Gamma_here = Gamma;
P.D   = create_demand(Gamma_here,Gamma_here*GammaH_ratio,0,mass);
tdlist{k} = compute_td(nfa0,T,P,shock,tdlist{k-1});k=k+1;

% learning: reduce Gamma over time
Gamma_ini = Gamma_here; Gamma_final = 0.1;
share_learned = 0.75; % how many have learned after 5 periods (.75 = no of intermediaries quadruples!)
learning = @(t) exp(log(1-share_learned)*t/5); % function between 1 and 0 that multiplies Gamma_ini and therefore can lower it over time!
learning_p = @(t) log(1-share_learned)/5*exp(log(1-share_learned)*t/5); % derivative
%     learning = @(t) Gamma_final/Gamma_ini + (Gamma_ini - Gamma_final)/Gamma_ini*(t<=5).*(1-t/5); % function between 1 and 0 that multiplies Gamma_ini and therefore can lower it over time!
%     learning_p = @(t) (Gamma_ini - Gamma_final)/Gamma_ini*(t<=5)*(-1/5); % derivative
P.D   = create_demand(Gamma_ini,Gamma_ini*GammaH_ratio,0,mass,learning,learning_p);
tdlist{k} = compute_td(nfa0,T,P,shock,tdlist{k-1});k=k+1;


%%% plot it
Z = Z_common;
Z.plot_name      = 'learning';
Z.colorlist      = {'black!70','green','blue'};
Z.stylelist      = {'solid','solid','solid'};
Z.legend_        = {'No intervention','No learning','Learning'};
Z.max_T          = 10;
Z.plot_name = [Z.plot_name suffix]; plot_this(Z,tdlist,H);







