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

Published with MATLAB® R2020a