Contents

function out = negative_paper()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code simulates the Keen model with monetary policy as
% specified in Grasselli and Lipton (2019)
%
% authors: M. Grasselli and A. Lipton
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
clc
warning off

Base Parameters - Table 2

a = -0.0401;    % constant parameter in phillips curve
b_phil = 0.0001;% coefficient phillips curve
c = -0.0065;    % constant parameter in investment function
d_i = -5;         % affine parameter in exponent of investment function
e = 20;         % coefficient in the exponent of investment function
kr= 0.1;        % capital ratio for banks (loan/equity)
markup = 1.3;   % markup factor
alpha = 0.025;  % growth rate in productivity
beta =0.02;     % growth rate in population
gamma=0.8;      % money illusion coefficient
delta=0.03;     % depreciation rate
etap=0.35;      % inflation relaxation parameter
nu=3;           % capital to output ratio

Other definitions

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

Functions of interest

inflation = @(x) etap*(markup*x-1);              % inflation, equation (10)

phillips = @(x) a+b_phil./(1-x).^2;                  % Phillips curve, equation (24)

investment= @(x) c+exp(d_i+e.*x);                 % investment function, equation (25)

Example 1: good equilibirum of Keen model without policy rate

tax=0;        % tax to income ratio
g=0;          % government spending to income ratio
etag=0;       % speed of adjustment of rho
etar=0;       % speed of adjustment of rg
deltar=0.03;  % interest rate spread

T=100;

omega0=0.8;
lambda0=0.9;
l0=0.6;
rho0=0;
rg0=0;
b0=0;
Y0=1;

y0=[convert([omega0,lambda0]),l0,rho0,rg0,b0];

[tK,yK] = ode15s(@(t,y) main(y), [0 T], y0, options);
sol = retrieve([yK(:,1),yK(:,2)]);
yK(:,1)=sol(:,1);
yK(:,2)=sol(:,2);
Y_output = Y0*yK(:,2)/lambda0.*exp((alpha+beta)*tK);


% Figure 2
figure
subplot(2,3,1)
plot(tK,yK(:,1))
xlabel('t')
ylabel('\omega')

subplot(2,3,2)
plot(tK,yK(:,2))
xlabel('t')
ylabel('\lambda')

subplot(2,3,3)
plot(tK,yK(:,3))
xlabel('t')
ylabel('ell')

subplot(2,3,4)
plot(tK,Y_output)
xlabel('t')
ylabel('Y')

inflation_sol=inflation(yK(:,1));
subplot(2,3,5)
plot(tK,inflation_sol)
xlabel('t')
ylabel('inflation')

subplot(2,3,6)
plot(yK(:,1),yK(:,2))
xlabel('\omega')
ylabel('\lambda')

Example 2: bad equilibrium of Keen model without policy rate

tax=0;        % tax to income ratio
g=0;          % government spending to income ratio
etag=0;       % speed of adjustment of rho
etar=0;       % speed of adjustment of rg
deltar=0.03;  % interest rate spread

T=100;

omega0=0.8;
lambda0=0.9;
l0=6;         % NOTE: this is incorrectly listed in the paper as l0 = 0.6
rho0=0;
rg0=0;
b0=0;
Y0=1;

y0=[convert([omega0,lambda0]),l0,rho0,rg0,b0];

[tK,yK] = ode15s(@(t,y) main(y), [0 T], y0, options);
sol = retrieve([yK(:,1),yK(:,2)]);
yK(:,1)=sol(:,1);
yK(:,2)=sol(:,2);
Y_output = Y0*yK(:,2)/lambda0.*exp((alpha+beta)*tK);

% Figure 3
figure
subplot(2,3,1)
plot(tK,yK(:,1))
xlabel('t')
ylabel('\omega')

subplot(2,3,2)
plot(tK,yK(:,2))
xlabel('t')
ylabel('\lambda')

subplot(2,3,3)
plot(tK,yK(:,3))
xlabel('t')
ylabel('ell')

subplot(2,3,4)
plot(tK,Y_output)
xlabel('t')
ylabel('Y')

inflation_sol=inflation(yK(:,1));
subplot(2,3,5)
plot(tK,inflation_sol)
xlabel('t')
ylabel('inflation')

subplot(2,3,6)
plot(yK(:,1),yK(:,2))
xlabel('\omega')
ylabel('\lambda')

Example 3: good equilibrium with policy rate

tax=0;        % tax to income ratio
g=0.2;          % government spending to income ratio
etag=0.2;       % speed of adjustment of rho
etar=0.1;       % speed of adjustment of rg
deltar=0.03;  % interest rate spread

T=100;

omega0=0.8;
lambda0=0.9;
l0=0.6;
rho0=0;
rg0=0;
b0=0.4;
Y0=1;

y0=[convert([omega0,lambda0]),l0,rho0,rg0,b0];

[tK,yK] = ode15s(@(t,y) main(y), [0 T], y0, options);
sol = retrieve([yK(:,1),yK(:,2)]);
yK(:,1)=sol(:,1);
yK(:,2)=sol(:,2);
Y_output = Y0*yK(:,2)/lambda0.*exp((alpha+beta)*tK);

% Figure 4
figure
subplot(3,3,1)
plot(tK,yK(:,1))
xlabel('t')
ylabel('\omega')

subplot(3,3,2)
plot(tK,yK(:,2))
xlabel('t')
ylabel('\lambda')

subplot(3,3,3)
plot(tK,yK(:,3))
xlabel('t')
ylabel('ell')

subplot(3,3,4)
plot(tK,yK(:,4))
xlabel('t')
ylabel('\rho')

subplot(3,3,5)
plot(tK,yK(:,5))
xlabel('t')
ylabel('rg')

subplot(3,3,6)
plot(tK,yK(:,6))
xlabel('t')
ylabel('b')

subplot(3,3,7)
plot(tK,Y_output)
xlabel('t')
ylabel('Y')

inflation_sol=inflation(yK(:,1));
subplot(3,3,8)
plot(tK,inflation_sol)
xlabel('t')
ylabel('inflation')

subplot(3,3,9)
plot(yK(:,1),yK(:,2))
xlabel('\omega')
ylabel('\lambda')

Example 4: stabilizaing policy rate with high debt levels

tax=0;        % tax to income ratio
g=0.2;          % government spending to income ratio
etag=0.2;       % speed of adjustment of rho
etar=0.1;       % speed of adjustment of rg
deltar=0.03;  % interest rate spread

T=100;

omega0=0.8;
lambda0=0.9;
l0=6;
rho0=0;
rg0=0;
b0=0.4;
Y0=1;

y0=[convert([omega0,lambda0]),l0,rho0,rg0,b0];

[tK,yK] = ode15s(@(t,y) main(y), [0 T], y0, options);
sol = retrieve([yK(:,1),yK(:,2)]);
yK(:,1)=sol(:,1);
yK(:,2)=sol(:,2);
Y_output = Y0*yK(:,2)/lambda0.*exp((alpha+beta)*tK);

% Figure 5
figure
subplot(3,3,1)
plot(tK,yK(:,1))
xlabel('t')
ylabel('\omega')

subplot(3,3,2)
plot(tK,yK(:,2))
xlabel('t')
ylabel('\lambda')

subplot(3,3,3)
plot(tK,yK(:,3))
xlabel('t')
ylabel('ell')

subplot(3,3,4)
plot(tK,yK(:,4))
xlabel('t')
ylabel('\rho')

subplot(3,3,5)
plot(tK,yK(:,5))
xlabel('t')
ylabel('rg')

subplot(3,3,6)
plot(tK,yK(:,6))
xlabel('t')
ylabel('b')

subplot(3,3,7)
plot(tK,Y_output)
xlabel('t')
ylabel('Y')

inflation_sol=inflation(yK(:,1));
subplot(3,3,8)
plot(tK,inflation_sol)
xlabel('t')
ylabel('inflation')

subplot(3,3,9)
plot(yK(:,1),yK(:,2))
xlabel('\omega')
ylabel('\lambda')

Example 5: stabilizaing policy rate with higher debt levels

tax=0;          % tax to income ratio
g=0.2;          % government spending to income ratio
etag=0.2;       % speed of adjustment of rho
etar=0.1;       % speed of adjustment of rg
deltar=0.03;    % interest rate spread

T=100;

omega0=0.8;
lambda0=0.9;
l0=8;
rho0=0;
rg0=0;
b0=0.4;
Y0=1;

y0=[convert([omega0,lambda0]),l0,rho0,rg0,b0];

[tK,yK] = ode15s(@(t,y) main(y), [0 T], y0, options);
sol = retrieve([yK(:,1),yK(:,2)]);
yK(:,1)=sol(:,1);
yK(:,2)=sol(:,2);
Y_output = Y0*yK(:,2)/lambda0.*exp((alpha+beta)*tK);

% Figure 6
figure
subplot(3,3,1)
plot(tK,yK(:,1))
xlabel('t')
ylabel('\omega')

subplot(3,3,2)
plot(tK,yK(:,2))
xlabel('t')
ylabel('\lambda')

subplot(3,3,3)
plot(tK,yK(:,3))
xlabel('t')
ylabel('ell')

subplot(3,3,4)
plot(tK,yK(:,4))
xlabel('t')
ylabel('\rho')

subplot(3,3,5)
plot(tK,yK(:,5))
xlabel('t')
ylabel('rg')

subplot(3,3,6)
plot(tK,yK(:,6))
xlabel('t')
ylabel('b')

subplot(3,3,7)
plot(tK,Y_output)
xlabel('t')
ylabel('Y')

inflation_sol=inflation(yK(:,1));
subplot(3,3,8)
plot(tK,inflation_sol)
xlabel('t')
ylabel('inflation')

subplot(3,3,9)
plot(yK(:,1),yK(:,2))
xlabel('\omega')
ylabel('\lambda')

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)

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

function new = convert(old)
    % This function converts
    % [omega, lambda]
    % to
    % [log_omega, tan_lambda]


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

    new(:,1) = log(old(:,1));
    new(:,2) = tan((old(:,2)-0.5)*pi);
end

function old = retrieve(new)
    % This function retrieves [omega, lambda]
    % from
    % [log_omega, tan_lambda]

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

    old(:,1) = exp(new(:,1));
    old(:,2) = atan(new(:,2))/pi+0.5;
end

function f = main(y)    % solves the main dynamical system
    f = zeros(6,1);

    log_omega = y(1);
    tan_lambda = y(2);
    l = y(3);
    rho= y(4);
    rg=y(5);
    b=y(6);

    lambda = atan(tan_lambda)/pi+0.5;
    omega = exp(log_omega);
    profit = 1-omega-tax-(rg+deltar)*l;

    growth = investment(profit)/nu-delta;    % real growth rate

    % Equation (13)
    f(1) = phillips(lambda)-(1-gamma)*inflation(omega)-alpha; %d(log_omega)/dt
    f(2) = (1+tan_lambda^2)*pi*lambda*(growth-alpha-beta); %d(tan_lambda)/dt
    f(3) = (rg+deltar-growth-inflation(omega))*l+omega+tax+investment(profit)-1; %d(l)/dt
    f(4) = etag*(growth-alpha-beta);     %d(rho)/dt
    f(5) = etar*(rho-rg);                %d(rg)/dt

    % Equation (18)
    f(6) = (g-tax)+b*(rg-inflation(omega)-growth);   %db/dt

end

end


Published with MATLAB® R2017a