## Contents

function out = ponzi()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code simulates the Keen model with Ponzi financing
% as specified in Grasselli and Costa Lima (2012)
%
%
% authors: B. Costa Lima and M. Grasselli
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
clc
warning off


Equations (24), (26), (63), (65) and (90)

nu = 3;              % capital-to-output ratio
alpha = 0.025;       % productivity growth rate
beta = 0.02;         % population growth rate
delta = 0.01;        % depreciation rate
phi0 = 0.04/(1-0.04^2); %Phillips curve parameter,  Phi(0.96)=0
phi1 = 0.04^3/(1-0.04^2); %Phillips curve parameter, Phi(0)=-0.04
r = 0.03;            %real interest rate
kappa0 = -0.0065 ;     % investment function parameter
kappa1 = exp(-5);      % investment function parameter
kappa2 = 20;           % investment function parameter
psi0 =-0.25;           % speculation function parameter
psi1 =0.25*exp(-0.36); % speculation function parameter
psi2 =12;


## Other definitions

TOL = 1E-7;
options = odeset('RelTol', TOL);
txt_format = '%3.3g';


## Functions Definition

fun_Phi = @(x) phi1./(1-x).^2- phi0;           % Phillips curve, equation (25)
fun_Phi_prime = @(x) 2*phi1./(1-x).^3;         % derivative of Phillips curve
fun_Phi_inv = @(x) 1 - (phi1/(x+phi0)).^(1/2); % inverse of Phillips curve

fun_kappa = @(x) kappa0+kappa1.*exp(kappa2.*x);         % Investment function, equation (64)
fun_kappa_prime = @(x) kappa1.*kappa2.*exp(kappa2.*x);  % Derivative of investment function
fun_kappa_inv = @(x) log((x-kappa0)./kappa1)./kappa2;   % Inverse of investment function

fun_Psi = @(x) psi0+psi1.*exp(psi2.*x);         % speculation function, equation (89)


## Sample paths of the Keen model with speculation

bar_pi_n_K = fun_kappa_inv(nu*(alpha+beta+delta));
bar_d_K = (fun_kappa(bar_pi_n_K)-bar_pi_n_K)/(alpha+beta);
bar_omega_K = 1-bar_pi_n_K-r*bar_d_K;
bar_lambda_K = fun_Phi_inv(alpha);

omega0 = 0.95;
lambda0 = 0.9;
d0 = 0;
p0 =0.1;
Y0 = 100;
T = 200;

[tK,yK] = ode15s(@(t,y) keen(y), [0 T], convert([omega0, lambda0, d0]), options);
yK = retrieve(yK);
Y_output = Y0*yK(:,2)/lambda0.*exp((alpha+beta)*tK);

[tP,yP] = ode15s(@(t,y) ponzi(y), [0 T], convert([omega0, lambda0, d0, p0]), options);
yP = retrieve(yP);
Y_output_ponzi = Y0*yP(:,2)/lambda0.*exp((alpha+beta)*tP);

% Figure 7
figure(1)
plot(yK(:,2), yK(:,1),'-k'), hold on;
plot(yP(:,2), yP(:,1),'--b')
xlabel('\lambda')
ylabel('\omega')
title(['\omega_0 = ', num2str(omega0, txt_format), ...
', \lambda_0 = ', num2str(lambda0, txt_format), ...
', d_0 = ', num2str(d0, txt_format), ...
', p_0 = ', num2str(p0, txt_format), ...
', Y_0 = ', num2str(Y0, txt_format)])
legend('No Speculation ', 'Ponzi Financing','Location','NorthWest')
print('keen_ponzi_phase','-depsc')

% Figure 8
figure(2)
fig = figure;
left_color = [0 0 1];
right_color = [0.5 0.7 0.2];
set(fig,'defaultAxesColorOrder',[left_color; right_color]);
title(['\omega_0 = ', num2str(omega0, txt_format), ...
', \lambda_0 = ', num2str(lambda0, txt_format), ...
', d_0 = ', num2str(d0, txt_format), ...
', p_0 = ', num2str(p0, txt_format), ...
', Y_0 = ', num2str(Y0, txt_format)])
subplot(3,1,1)
yyaxis right
plot(tP, yP(:,3),'--','color',[0.5 0.7 0.2]')
ylabel('d with Ponzi')

yyaxis left
plot(tK, yK(:,3),'-b')
ylabel('d')
legend('No Speculation ', 'Ponzi Financing','Location','NorthEast')

subplot(3,1,2)
yyaxis right
plot(tP, Y_output_ponzi,'--','color',[0.5 0.7 0.2]')
ylabel('Y with Ponzi')

yyaxis left
plot(tK, Y_output,'-b')
ylabel('Y')
legend('No Speculation ', 'Ponzi Financing','Location','NorthEast')

subplot(3,1,3)

plot(tP, yP(:,4),'-b')
ylabel('p')
print('keen_ponzi_time','-depsc')


## Auxiliary functions

Instead of integrating the system in terms of omega, lambda and d, we decided to use:

log_omega = log(omega), tan_lambda = tan((lambda-0.5)*pi), and pi_n=1-omega-r*d log_p = log(p)

Experiments show that the numerical methods work better with these variables.

function new = convert(old)
% This function converts
% [omega, lambda, d, p]
% to
% [log_omega, tan_lambda, pi_n, log_p]
% It also work with the first 2 or 3 elements alone.
% each of these variables might also be a time-dependent vector

n = size(old,2); %number of variables
new = zeros(size(old));

new(:,1) = log(old(:,1));
new(:,2) = tan((old(:,2)-0.5)*pi);
if n>2
new(:,3) = 1-old(:,1)-r*old(:,3);
if n==4
new(:,4) = log(old(:,4));
end
end
end

function old = retrieve(new)
% This function retrieves [omega, lambda, d, p]
% from
% [log_omega, tan_lambda, pi_n, log_p]
%
%  It also work with the first 2 or 3 elements alone.
% each of these variables might also be a time-dependent vector

n = size(new,2); %number of variables
old = zeros(size(new));

old(:,1) = exp(new(:,1));
old(:,2) = atan(new(:,2))/pi+0.5;
if n>2
old(:,3) = (1-old(:,1)-new(:,3))/r;
if n==4
old(:,4) = exp(new(:,4));
end
end
end

function f = keen(y)
f = zeros(3,1);
log_omega = y(1);
tan_lambda = y(2);
pi_n = y(3);
lambda = atan(tan_lambda)/pi+0.5;
omega = exp(log_omega);
g_Y = fun_kappa(pi_n)/nu-delta;

% Equation (47)
f(1) = fun_Phi(lambda)-alpha;                                   %d(log_omega)/dt
f(2) = (1+tan_lambda^2)*pi*lambda*(g_Y-alpha-beta);             %d(tan_lambda)/dt
f(3) = -omega*f(1)-r*(fun_kappa(pi_n)-pi_n)+(1-omega-pi_n)*g_Y; %d(pi_n)/dt
end

function f = ponzi(y)
f = zeros(4,1);
log_omega = y(1);
tan_lambda = y(2);
pi_n = y(3);
log_p = y(4);
lambda = atan(tan_lambda)/pi+0.5;
omega = exp(log_omega);
p = exp(log_p);
g_Y = fun_kappa(pi_n)/nu-delta;

% Equation (74)
f(1) = fun_Phi(lambda)-alpha;                                       %d(log_omega)/dt
f(2) = (1+tan_lambda^2)*pi*lambda*(g_Y-alpha-beta);                 %d(tan_lambda)/dt
f(3) = -omega*f(1)-r*(fun_kappa(pi_n)-pi_n)+(1-omega-pi_n)*g_Y-r*p; %d(pi_n)/dt
f(4) = fun_Psi(g_Y) - g_Y;                                          %dp/dt
end

end