## Contents

function out = keen_1995()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code simulates the Keen model as specified in
% S. Keen. Finance and economic breakdown: Modeling Minsky’s “Financial Instability Hypothesis”.
% Journal of Post Keynesian Economics, 17(4):607–635, 1995.
%
% author: M. Grasselli
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
clc
warning off


## Parameters attribution and function definitions

% General parameters, Page 615, footnote 2
par_alpha = 0.015;       % productivity growth rate
par_beta = 0.035;        % population growth rate
par_gamma = 0.02;        % depreciation rate
par_nu = 3;              % capital-to-output ratio

% Phillips curve, page 615, Phi(lambda) = A/(1-lambda)^2-D
par_A = 0.0000641;         % Phi(0)=-0.04
par_D = 0.0400641;         % Phi(0.96)=0
fun_Phi = @(x) par_A./(1-x).^2- par_D;           % Phillips curve
fun_Phi_inv = @(y) 1 - sqrt(par_A./(y+par_D));   % inverse of Phillips curve

% Investment function, page 616, kappa(pi/nu) = E/(F-G*pi/nu)^2 - H
par_E = 0.0175 ;
par_F = 0.53;
par_G = 6;
par_H = 0.065;
% Investment function should have been capped at 1 to prevent consumption going negative
fun_kappa = @(x) min(1,par_E./(par_F-(par_G/par_nu).*x).^2-par_H);  % investment function
fun_kappa_inv = @(y) (par_nu/par_G)*(par_F-sqrt(par_E/(y+par_H)));  % inverse of investment function

% Government subsidy function, page 625, footnote 4, G(1-lambda) = i/(j-k*(1-lambda))^2 - l
par_i = 0.05;
par_j = 1.2;
par_k = 4;
par_l = 0.05;
fun_Gamma = @(x) par_i./(par_j-par_k.*x).^2- par_l;

% Taxation function, page 625, footnote 4, T(pi) = m/(n-o*pi)^2 - p
par_m = 0.0175;
par_n = 0.83;
par_o = 1.6014;     % this is o = 5 in the paper, but does not match Figure 11
par_p = 0.039;
fun_Theta = @(x) par_m./(par_n-par_o.*x).^2- par_p;

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


## Basic Goodwin limit cycle, page 617

% Figure 1 in the paper
figure(1)

tiledlayout(1,2)
sgtitle('Behavioral functions for workers and capitalists')
nexttile
lambda_var=[0.9:0.001:0.99];
y1=fun_Phi(lambda_var);
plot(lambda_var,y1)
xlabel('Employment Rate')
ylabel('Rate of change of real wage')

nexttile
pi_var=[-0.1:0.001:0.25];
y2=fun_kappa(pi_var);
plot(pi_var,y2)
xlabel('Profit/Output')
ylabel('Investment/Output')

omega0 = 0.96;                                  % initial wage share, footnote 2
lambda0 = 0.90;                                 % initial employment rate, footnote 2
T = 100;                                        % time in years
lambda_bar = fun_Phi_inv(par_alpha);            % equilibrium employment rate
omega_bar =  1-par_nu*(par_alpha+par_beta+par_gamma);        % equilibrium wage share in original Goodwin model
omega_mod = 1-fun_kappa_inv(par_nu*(par_alpha+par_beta+par_gamma)); % equilibrium wage share in the modified Goodwin model with investment function

par1 = [par_nu,par_alpha,par_beta,par_gamma,par_A,par_D];
par2 = [par_nu,par_alpha,par_beta,par_gamma,par_A,par_D,par_E,par_F,par_G,par_H];

% This uses the function 'goodwin' defined below to find the trajectories of the Goodwin model
[tG,yG] = ode15s(@(t,y) goodwin(y,par1), [0 T], [omega0, lambda0], options);
[tGmod,yGmod] = ode15s(@(t,y) goodwin_mod(y,par2), [0 T], [omega0, lambda0], options);

% Figure 2 in the paper
figure(2)
plot(tG, yG(:,2), 'k', tG, yG(:,1), 'r')
xlabel('time')
%ylabel('\omega,\lambda')
legend('Employment','Wage Share')
title(['Wage share and employment, basic Goodwin model'])

% Figure 3 in the paper
figure(3)
plot(yG(:,1), yG(:,2)), hold on;
plot(omega_bar, lambda_bar, 'ro')
plot(yGmod(:,1), yGmod(:,2))
plot(omega_mod, lambda_bar, 'rx')
legend('Basic Goodwin Cycle','Basic Goodwin Equilibrium','Investment Function Cycle',...
'Investment Function Equilibrium')
xlabel('Wage Share')
ylabel('Employment')
title(['Cyclical and equilibrium time paths'])


## Base rate variations, page 620

% Low base rate

par_xi = 0.036;         % base value for real interest rate (the paper says it should be below 0.046), this
% was reversed engineered to match the equilibrium
% value for b in Figure 5

par_phi = 0;

% Numerical Values of Interest
kappa_eq = par_nu*(par_alpha+par_beta+par_gamma);                % investment at interior equilibrium
pi_eq=fun_kappa_inv(par_nu*(par_alpha+par_beta+par_gamma));      % profit at interior equilibrium
d_eq = (fun_kappa(pi_eq)-pi_eq-par_nu*par_gamma)/(par_alpha+par_beta-par_xi);   % debt ratio at interior equilibirum
omega_eq = 1-pi_eq-par_xi*d_eq;                                  % wage share at interior equilibrium
lambda_eq=fun_Phi_inv(par_alpha);                                % employment at interior equilibirum
b_eq = par_xi*d_eq;

omega0 = 0.96;                                 % initial wage share, footnote 2
lambda0 = 0.90;                                % initial employment rate, footnote 2
d0 = 0.0;                                      % initial debt, remark about zero banker's share on page 620

y0 = [convert([omega0, lambda0]),d0];          % change of variable for first to components

T = 200;

par3 = [par_nu,par_alpha,par_beta,par_gamma,par_xi,par_phi,par_A,par_D,par_E,par_F,par_G,par_H];

% This uses the function 'keen_wrong' defined below to find the trajectories for the Keen
% model as specified on page 616
[tK,yK] = ode15s(@(t,y) keen_wrong(y,par3), [0 T], y0, options);
yKnew = retrieve([yK(:,1),yK(:,2)]);
yK(:,1) = yKnew(:,1);
yK(:,2) = yKnew(:,2);

% Figure 4 in the paper
figure(4)
plot(tK, yK(:,2),tK,yK(:,1))
xlabel('time')
legend('Employment','Wage Share')
title(['Employment and wage share at low interest'])

% Figure 5 in the paper
figure(5)
profit = 1-yK(:,1)-par_xi*yK(:,3);
invest = fun_kappa(profit);
bank = par_xi*yK(:,3);
plot(tK,invest,tK,profit,tK,bank)
legend('Investment','Profit Share','Bank Share')
ylim([0 0.4])
title(['Profit, investment, and bank share at low interest'])

% Figure 6 in the paper
figure(6)
plot3(yK(:,2),yK(:,1),bank)
xlabel('Employment')
ylabel('Wages')
zlabel('Bank Share')
title(['Income and interest stability at low interest'])

% High base rate

par_xi = 0.04626;       % base value for real interest rate (the paper says it should be above 0.046), this
% was reversed engineered to match the dynamics for
% omega and lambda on Figure 7

par_phi = 0;

% Numerical Values of Interest
kappa_eq = par_nu*(par_alpha+par_beta+par_gamma);                % investment at interior equilibrium
pi_eq=fun_kappa_inv(par_nu*(par_alpha+par_beta+par_gamma));      % profit at interior equilibrium
d_eq = (fun_kappa(pi_eq)-pi_eq-par_nu*par_gamma)/(par_alpha+par_beta-par_xi);   % debt ratio at interior equilibirum
omega_eq = 1-pi_eq-par_xi*d_eq;                                  % wage share at interior equilibrium
lambda_eq=fun_Phi_inv(par_alpha);                                % employment at interior equilibirum
b_eq = par_xi*d_eq;

omega0 = 0.96;                                 % initial wage share, footnote 2
lambda0 = 0.90;                                % initial employment rate, footnote 2
d0 = 0.0;                                      % page 620

y0 = [convert([omega0, lambda0]),d0];

T = 1000;

par4 = [par_nu,par_alpha,par_beta,par_gamma,par_xi,par_phi,par_A,par_D,par_E,par_F,par_G,par_H];

[tK,yK] = ode15s(@(t,y) keen_wrong(y,par4), [0 T], y0, options);
yKnew = retrieve([yK(:,1),yK(:,2)]);
yK(:,1) = yKnew(:,1);
yK(:,2) = yKnew(:,2);

% Figure 7 in the paper
figure(7)

tiledlayout(2,1)
nexttile
plot(tK, yK(:,2),tK,yK(:,1))
xlabel('time')
legend('Employment','Wage Share')
title(['Instability at "high" interest'])

% Figure not shown in the paper, but mentioned on page 622
nexttile
profit = 1-yK(:,1)-par_xi*yK(:,3);
invest = fun_kappa(profit);
bank = par_xi*yK(:,3);
plot(tK,invest,tK,profit,tK,bank)
legend('Investment','Profit Share','Bank Share')
ylim([0 0.4])
title(['Profit, investment, and bank share at high interest'])

% Figure 8 in the paper
figure(8)
plot3(yK(:,2),yK(:,1),bank)
xlabel('Employment')
ylabel('Wage Share')
zlabel('Bank Share')
title(['Instability at "high" interest'])


## Debt sensitivity variations, page 623

% Low base rate, high sensitivity

par_xi = 0.0459;        % base value for real interest rate (the paper says it should be below 0.046), this
% was reversed engineered to match the dynamics for
% omega and lambda in Figure 9

par_phi = 0.0022;       % tried different values to match the dynamics in Figures 9 and 10

omega0 = 0.96;                                 % initial wage share, footnote 2
lambda0 = 0.90;                                % initial employment rate, footnote 2
d0 = 0.0;                                      % page 620

y0 = [convert([omega0, lambda0]),d0];

T = 300;

par5 = [par_nu,par_alpha,par_beta,par_gamma,par_xi,par_phi,par_A,par_D,par_E,par_F,par_G,par_H];

[tK,yK] = ode15s(@(t,y) keen_wrong(y,par5), [0 T], y0, options);
yKnew = retrieve([yK(:,1),yK(:,2)]);
yK(:,1) = yKnew(:,1);
yK(:,2) = yKnew(:,2);

% Figure 9 in the paper
figure(9)
tiledlayout(2,1)
nexttile
plot(tK, yK(:,2),tK,yK(:,1))
xlabel('time')
legend('Employment','Wage Share')
title(['Debt sensitivity driven crisis'])

% Not shown in the paper
nexttile
profit = 1-yK(:,1)-(par_xi+par_phi*yK(:,3)).*yK(:,3);
invest = fun_kappa(profit);
bank = (par_xi+par_phi*yK(:,3)).*yK(:,3);
plot(tK,invest,tK,profit,tK,bank)
legend('investment','profit','bank share')
ylim([0 0.4])
title(['Profit, investment, and bank share at low interest, high sensitivity'])

% Figure 10 in the paper
figure(10)
plot3(yK(:,2),yK(:,1),bank)
xlabel('Emplyment')
ylabel('Wage Share')
zlabel('Bank Share')
title(['Debt sensitivity crisis'])


## Minskian government, page 625

% Figure 11 in the paper
figure(11)

lambda_var = [0.8:0.01:1];
y3 = fun_Gamma(1-lambda_var);

pi_var = [-0.3:0.01:0.3];
y4 = fun_Theta(pi_var);

tiledlayout(1,2)
sgtitle('Government behavioral functions')
nexttile
plot(lambda_var,y3)
xlabel('Employment rate')
ylabel('Rate of change of government subsidy')

nexttile
plot(pi_var,y4)
xlabel('Gross Profit/Output')
ylabel('Rate of change of taxes')

omega0 = 0.96;                                 % initial wage share, footnote 2
lambda0 = 0.90;                                % initial employment rate, footnote 2
dk0 = 0;
dg0 = 0;
g0 = 0;
t0 = 0;

par_xi = 0.05;              % value provided in Figure 12
par_phi = 0.005;            % value provided in Figure 12

y0 = [convert([omega0, lambda0]),dk0,dg0,g0,t0];

T = 350;

par6 = [par_nu,par_alpha,par_beta,par_gamma,par_xi,par_phi,par_A,par_D,par_E,par_F,par_G,par_H,...
par_i,par_j,par_k,par_l,par_m,par_n,par_o,par_p];

% This uses the function 'keen_government' defined below to find the trajectories of the
% Keen model with government as specified on page 628
[tG,yG] = ode15s(@(t,y) keen_government(y,par6), [0 T], y0, options);
yGnew = retrieve([yG(:,1),yG(:,2)]);
yG(:,1) = yGnew(:,1);
yG(:,2) = yGnew(:,2);

% Figure 12 in the paper
figure(12)
tiledlayout(3,1)
nexttile
plot(tG, yG(:,2),tG,yG(:,1))
xlabel('time')
legend('Employment','Wage Share')
title(['Persistent cycles with government'])

% Not shown in the paper
nexttile
profit = 1-yG(:,1)-(par_xi+par_phi*yG(:,3)).*yG(:,3)+yG(:,5)-yG(:,6);
gross_profit = 1-yG(:,1);
invest = fun_kappa(profit);
bank = (par_xi+par_phi*yG(:,3)).*yG(:,3);
plot(tG,gross_profit,tG,invest,tG,profit,tG,bank)
legend('gross profit','investment','net profit','bank share')
title(['Profit, investment, and bank share with government'])
ylim([-1,1])

% Not shown in the paper
nexttile
plot(tG, yG(:,4),tG,yG(:,5),tG,yG(:,6))
legend('Government debt','subsidy','tax')
title(['Government debt, subsidy, and tax ratios'])
ylim([-1,1])

% Figure 13 in the paper
figure(13)
tiledlayout(2,1)
sgtitle(['The stabilizing effect of countercyclical government'])
nexttile
plot(yG(:,2), yG(:,1))
ylabel('Wage Share')
xlabel('Employment')
title(['Private Sector'])

nexttile
net_government = yG(:,5)-yG(:,6);
plot(yG(:,2), net_government)
ylabel('Net Government Spending Share')
xlabel('Employment')
title(['Public Sector'])

% Figure 14 in the paper
figure(14)
plot3(yG(:,2),yG(:,1),bank)
xlabel('Employment')
ylabel('Wages')
zlabel('Bank Share')
title(['Complexity with government, but no breakdown'])

% Figure 15 in the paper
figure(15)
step=round(length(tG)/7);    % This corresponds to 50 periods
tiledlayout(3,2)
sgtitle(['Cyclical stability with government'])
nexttile
plot(yG(1:2*step,2),yG(1:2*step,1))
xlabel('Employment')
ylabel('Wages')

nexttile
plot(yG(2*step:3*step,2),yG(2*step:3*step,1))
xlabel('Employment')
ylabel('Wages')

nexttile
plot(yG(3*step:4*step,2),yG(3*step:4*step,1))
xlabel('Employment')
ylabel('Wages')

nexttile
plot(yG(4*step:5*step,2),yG(4*step:5*step,1))
xlabel('Employment')
ylabel('Wages')

nexttile
plot(yG(5*step:6*step,2),yG(5*step:6*step,1))
xlabel('Employment')
ylabel('Wages')

nexttile
plot(yG(6*step:7*step,2),yG(6*step:7*step,1))
xlabel('Employment')
ylabel('Wages')


## Auxiliary functions

% Original Goodwin model
function f = goodwin(y,par)
f = zeros(2,1);
omega = y(1);
lambda = y(2);

nu = par(1);
alpha = par(2);
beta = par(3);
gamma = par(4);

A = par(5);
D = par(6);

phillips = A./(1-lambda)^2- D;

f(1) = omega*(phillips-alpha);                        % d(omega)/dt
f(2) = lambda*((1-omega)/nu-alpha-beta-gamma);        % d(lambda)/dt
end

% Goodwin model with nonlinear investment function
function f = goodwin_mod(y,par)
f = zeros(2,1);
omega = y(1);
lambda = y(2);

nu = par(1);
alpha = par(2);
beta = par(3);
gamma = par(4);

A = par(5);
D = par(6);

E = par(7);
F = par(8);
G = par(9);
H = par(10);

phillips = A./(1-lambda)^2- D;
kappa = E./(F-(G/nu).*(1-omega)).^2-H;

f(1) = omega*(phillips-alpha);                    % d(omega)/dt
f(2) = lambda*(kappa/nu-alpha-beta-gamma);        % d(lambda)/dt
end

% For the Keen model, 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)    (note pi = 3.15 here, NOT profits)

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

function new = convert(old,r)
% 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,r)
% 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

% Keen model as specified on page 617
% Notice that this is wrong for two reasons.
%
% First, the investment defined on page 615 should be gross, rather than net investment.
% As a result, there is an extra term nu*gamma in the differential equation for d.

% Second, there should be no term (rD) in the equation for dD/dt on page 616, as this
% is already included in Pi = 1-W-rD. As a result, there is an extra term r*d in the differential equation for d.

% Both mistakes are corrected in later versions of the model.

function f = keen_wrong(y,par)
f = zeros(3,1);
log_omega = y(1);
tan_lambda = y(2);
d = y(3);
lambda = atan(tan_lambda)/pi+0.5;   % pi = 3.14 here
omega = exp(log_omega);

nu = par(1);
alpha = par(2);
beta = par(3);
gamma = par(4);
xi = par(5);
phi = par(6);

A = par(7);
D = par(8);

E = par(9);
F = par(10);
G = par(11);
H = par(12);

r = xi+phi*d;

pi_n = 1-omega-r*d;                      % we use pi_n here to avoid confusion with pi = 3.14
phillips = A./(1-lambda)^2- D;
kappa = E./(F-(G/nu).*pi_n).^2-H;

g_Y = kappa/nu-gamma;

f(1) = phillips-alpha;                                     %d(log_omega)/dt
f(2) = (1+tan_lambda^2)*pi*lambda*(g_Y-alpha-beta);        %d(tan_lambda)/dt (notice here pi = 3.14)
f(3) = r*d-pi_n+(nu-d)*(kappa/nu-gamma);                   %dd/dt

end

% Keen model with government as specified on page 628

% Notice that the equation for dD/dt on page 627 is correct (as opposed
% to the incorrect equation for dD/dt on page 616), because in this case
% Pi = 1-W does not include any interest payment.

function f = keen_government(y,par)
f = zeros(6,1);
log_omega = y(1);
tan_lambda = y(2);
dk = y(3);
dg = y(4);
g = y(5);
t = y(6);

lambda = atan(tan_lambda)/pi+0.5;   % pi = 3.14 here
omega = exp(log_omega);

nu = par(1);
alpha = par(2);
beta = par(3);
gamma = par(4);
xi = par(5);
phi = par(6);

A = par(7);
D = par(8);

E = par(9);
F = par(10);
G = par(11);
H = par(12);

i = par(13);
j = par(14);
k = par(15);
l = par(16);

m = par(17);
n = par(18);
o = par(19);
p = par(20);

r = xi+phi*(dk+dg);

pi_n = 1-omega-r.*dk+g-t;
phillips = A./(1-lambda)^2- D;
kappa = E./(F-(G/nu).*pi_n).^2-H;
Gamma = i./(j-k.*(1-lambda)).^2-l;
Theta = m./(n-o.*pi_n).^2-p;                % using pi_n as the argument for the tax function gives the correct simulations
%  Theta = m./(n-o.*(1-omega)).^2-p;           % using pi = 1-omega (as it is claimed in the paper) does not work

g_Y = kappa/nu-gamma;

f(1) = phillips-alpha;                                     % d(log_omega)/dt
f(2) = (1+tan_lambda^2)*pi*lambda*(g_Y-alpha-beta);        % d(tan_lambda)/dt
f(3) = r*dk-(1-omega)+t-g+(nu-dk)*(kappa/nu-gamma);        % dd_k/dt
f(4) = r*dg-dg*(kappa/nu-gamma)+g-t;                       % dd_g/dt
f(5) = Gamma - g*(kappa/nu-gamma);                         % dg/dt
f(6) = Theta - t*(kappa/nu-gamma);                         % d(t)/dt
end

end