## Contents

function out = keen()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code simulates the Keen model as specified in Grasselli and Costa Lima (2012)
%
%
% authors: B. Costa Lima and M. Grasselli
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
clc
warning off

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

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

## 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

## Numerical Values of Interest

kappa_eq = nu*(alpha+beta+delta)                % investment at interior equilibrium
pi_eq=fun_kappa_inv(nu*(alpha+beta+delta))      % profit at interior equilibrium, equation (41)
d_eq = (fun_kappa(pi_eq)-pi_eq)/(alpha+beta)    % debt ratio at interior equilibirum, equation (42)
omega_eq = 1-pi_eq-r*d_eq                       % wage share at interior equilibrium, equation (42)
lambda_eq=fun_Phi_inv(alpha)                    % employment at interior equilibirum, equation (42)

d=[];
syms d
d_sols(1) = vpasolve(d*(r+delta-fun_kappa(1-r*d)/nu)+fun_kappa(1-r*d)-1,d) ;        % solution to equation (40) given in equation (66)
d_sols(2) = vpasolve(d*(r+delta-fun_kappa(1-r*d)/nu)+fun_kappa(1-r*d)-1,[20 40]) ;  % solution to equation (40) given in equation (66)
d_sols
J = @(x,y,z) [fun_Phi(y)-alpha, x.*fun_Phi_prime(y), 0;
-y.*fun_kappa_prime(1-x-r*z)/nu, (fun_kappa(1-x-r*z)-nu*(alpha+beta+delta))/nu, -r*y.*fun_kappa_prime(1-x-r*z)/nu;
((z-nu).*fun_kappa_prime(1-x-r*z)+nu)/nu, 0, (nu*(r+delta)-fun_kappa(1-x-r*z)+r*(z-nu).*fun_kappa_prime(1-x-r*z))/nu];

E1=eig(J(omega_eq,lambda_eq,d_eq)) % eigenvalues of the Jacobian at the interior equilibrium
E2=eig(J(0,0,d_sols(1)))           % eigenvalues of the Jacobian at the first d0 equilibrium
E3=eig(J(0,0,d_sols(2)))           % eigenvalues of the Jacobian at the second d0 equilibrium

cond59= fun_kappa_prime(pi_eq)/nu*(pi_eq-nu*delta)-(alpha+beta) % expression in brackets in condition (59)
kappa_eq =

0.1650

pi_eq =

0.1618

d_eq =

0.0702

omega_eq =

0.8361

lambda_eq =

0.9686

d_sols =

[ 2.999995114426015986641604626889, 30.767016452389481647975528982775]

E1 =

-0.0352 + 1.9581i
-0.0352 - 1.9581i
-0.0450 + 0.0000i

E2 =

-180122.65958170398076808544264108
-13/200
180122.11658158037574626271733576

E3 =

-13/200
-0.046692277832906022704348899916326
0.20619779407835991738362093529345

cond59 =

0.1057

## Sample paths of the Keen model

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.75;
lambda0 = 0.75;
d0 = 0.1;
Y0 = 100;
T = 300;

[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);

% Figure 3
figure(1)
plot(yK(:,1), yK(:,2)), hold on;
plot(bar_omega_K, bar_lambda_K, 'ro')
xlabel('\omega')
ylabel('\lambda')
title(['\omega_0 = ', num2str(omega0, txt_format), ...
', \lambda_0 = ', num2str(lambda0, txt_format), ...
', d_0 = ', num2str(d0, txt_format), ...
', Y_0 = ', num2str(Y0, txt_format)])
print('keen_phase_convergent','-depsc')

% Figure 4
figure(2)
yyaxis right
plot(tK, yK(:,1),tK,yK(:,2),tK,yK(:,3))
ylabel('\omega,\lambda,d')

yyaxis left
plot(tK, Y_output)
ylabel('Y')

title(['\omega_0 = ', num2str(omega0,txt_format), ...
', \lambda_0 = ', num2str(lambda0, txt_format), ...
', d_0 = ', num2str(d0, txt_format), ...
', Y_0 = ', num2str(Y0, txt_format)])
legend('Y','\omega', '\lambda','d')
print('keen_time_convergent','-depsc')

omega0 = 0.75;
lambda0 = 0.70;
d0 = 0.1;
Y0 = 100;
T = 300;

[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);

% Figure 5
figure(3)
yyaxis right
plot(tK, yK(:,1),tK,yK(:,2))
ylabel('\omega,\lambda')

yyaxis left
plot(tK,yK(:,3),tK, Y_output)
ylabel('d,Y')

title(['\omega_0 = ', num2str(omega0,txt_format), ...
', \lambda_0 = ', num2str(lambda0, txt_format), ...
', d_0 = ', num2str(d0, txt_format), ...
', Y_0 = ', num2str(Y0, txt_format)])
legend('d', 'Y','\omega','\lambda')
print('keen_time_divergent','-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
end