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
Parameters attribution
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