dual_keen

Base Parameters - Table 3
Functions of interest
Other definitions
Example 1
Example 2
Example 3
Auxiliary functions

function out = dual_keen()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code simulates the dual Keen model as specified in
% Giraud and Grasselli (2019)
%
%
% authors: G. Giraud and M. Grasselli
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc
warning off

Base Parameters - Table 3

nu = 3;          %capital to output ratio
alpha = 0.025;   %productivity growth rate
beta = 0.02;     %population growth rate
delta = 0.05;    %depreciation rate
r = 0.03;        %real interest rate
eta_p = 0.35;    % adjustment speed for prices
markup = 1.6;    % mark-up value
gamma = 0.8;     % (1-gamma) is the degree of money ilusion

phi0 = 0.04/(1-0.04^2); %Phillips curve parameter
phi1 = 0.04^3/(1-0.04^2); %Phillips curve parameter
c_minus= 0.03; % consumption share hard-lower bound
A_c = 0; %Consumption share asymptotic lower bound
B_c = 5; % Consumption share growth rate
C_c = 1;
K_c = 0.9; %Consumption share upper bound
nu_c = 0.2;
Q_c= ((K_c-A_c)/(c_minus-A_c))^nu_c-C_c;

Functions of interest

fun_consumption =@(x) max(c_minus,A_c + ((K_c-A_c)./(C_c+Q_c*exp(-B_c*x)).^(1/nu_c))); %----first specification : generalised logistic consumption function------%
fun_consumption_inv = @(x) -log((((K_c-A_c)./(x-A_c)).^nu_c-C_c)./Q_c)./B_c; % inverse of the consumption function over positive income

fun_Phi = @(x) phi1./(1-x).^2- phi0;
fun_Phi_prime = @(x) 2*phi1./(1-x).^3;
fun_Phi_inv = @(x) 1 - (phi1./(x+phi0)).^(1/2);

Other definitions

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

Example 1

% Parameters as in Table 3
nu = 3;          %capital to output ratio
alpha = 0.025;   %productivity growth rate
beta = 0.02;     %population growth rate
delta = 0.05;    %depreciation rate
r = 0.03;        %real interest rate
eta_p = 0.35;    % adjustment speed for prices
markup = 1.6;    % mark-up value
gamma = 0.8;     % (1-gamma) is the degree of money ilusion

% Numerical Values of Interest
c_eq=1-nu*(alpha+beta+delta)
eta_1 = fun_consumption_inv(c_eq)

A_w = eta_p*markup; % Coefficients in the quadratic equation A_w omega^2 + B_w omega + C_w = 0
B_w = alpha+beta-eta_p*(1+eta_1*markup);
C_w = eta_1*(eta_p-alpha-beta)+r*(eta_1-c_eq);

x=[];
syms x
omega_1 = vpasolve(A_w*x^2+B_w*x+C_w ==0,x)
i_omega_1 = eta_p*(markup*omega_1-1)
g_nominal_1 = alpha+beta+i_omega_1
lambda_1 = fun_Phi_inv(alpha+(1-gamma)*eta_p*(markup*omega_1-1))
d_1 = (omega_1-eta_1)./r

eta_1_calc=omega_1-r.*d_1; % verification step
c_eq_calc=fun_consumption(eta_1_calc); %verification step

% Sample paths of the dual Keen model

omega0 = 0.75;
lambda0 = 0.9;
d0 = 0.5;
Y0 = 100;
T = 200;

[tK,yK] = ode15s(@(t,y) dual_keen(y), [0 T], convert([omega0, lambda0, d0]), options);
yK = retrieve(yK);
Y_output = Y0*yK(:,2)/lambda0.*exp((alpha+beta)*tK);

n=length(yK(:,3))

omega_bar = yK(n,1)
lambda_bar = yK(n,2)
d_bar = yK(n,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('dual_keen_phase_convergent','-depsc')

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('dual_keen_time_convergent','-depsc')

c_eq =

    0.7150


eta_1 =

    0.6059

 
omega_1 =
 
 0.49291938520463741461070997693151
 0.65763195476319826831120345628563
 
 
i_omega_1 =
 
 -0.073965144285403047818002412918353
  0.018273894667391030254273935519954
 
 
g_nominal_1 =
 
 -0.028965144285403047818002412918353
  0.063273894667391030254273935519954
 
 
lambda_1 =
 
 0.96429092337337114586046395893964
 0.96945784620030614474363301296987
 
 
d_1 =
 
 -3.7663032540113713771463174119765
  1.7241157312739904128701318998274
 

n =

        2556


omega_bar =

    0.6567


lambda_bar =

    0.9696


d_bar =

    1.7244

Example 2

% Parameters as in Table 3 except for eta_p and gamma

nu = 3;         %capital to output ratio
alpha = 0.025;  %productivity growth rate
beta = 0.02;    %population growth rate
delta = 0.05;   %depreciation rate
r = 0.03;       %real interest rate
eta_p = 0.45;   % adjustment speed for prices - MODIFIED
markup = 1.6;   % mark-up value
gamma = .96;    % (1-gamma) is the degree of money ilusion - MODIFIED

% Numerical Values of Interest
c_eq=1-nu*(alpha+beta+delta)
eta_1 = fun_consumption_inv(c_eq)

A_w = eta_p*markup; % Coefficients in the quadratic equation A_w omega^2 + B_w omega + C_w = 0
B_w = alpha+beta-eta_p*(1+eta_1*markup);
C_w = eta_1*(eta_p-alpha-beta)+r*(eta_1-c_eq);

x=[];
syms x
omega_1 = vpasolve(A_w*x^2+B_w*x+C_w ==0,x)
i_omega_1 = eta_p*(markup*omega_1-1)
g_nominal_1 = alpha+beta+i_omega_1
lambda_1 = fun_Phi_inv(alpha+(1-gamma)*eta_p*(markup*omega_1-1))
d_1 = (omega_1-eta_1)./r

eta_1_calc=omega_1-r.*d_1; % verification step
c_eq_calc=fun_consumption(eta_1_calc); %verification step


% Sample paths of the dual Keen model

omega0 = 0.75;
lambda0 = 0.9;
d0 = 0.5;
Y0 = 100;
T = 200;

[tK,yK] = ode15s(@(t,y) dual_keen(y), [0 T], convert([omega0, lambda0, d0]), options);
yK = retrieve(yK);
Y_output = Y0*yK(:,2)/lambda0.*exp((alpha+beta)*tK);

n=length(yK(:,3))

omega_bar = yK(n,1)
lambda_bar = yK(n,2)
d_bar = yK(n,3)

figure(3)
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('dual_keen_phase_cycle','-depsc')

figure(4)
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('dual_keen_time_cycle','-depsc')
%

c_eq =

    0.7150


eta_1 =

    0.6059

 
omega_1 =
 
 0.51337660572143848603532778591854
 0.65503187710354006988977171337227
 
 
i_omega_1 =
 
 -0.08036884388056429005456399413865
 0.021622951514548850320635633628032
 
 
g_nominal_1 =
 
 -0.03536884388056429005456399413865
 0.066622951514548850320635633628032
 
 
lambda_1 =
 
 0.96780635616165993232539783359532
 0.96881832873512291925191143457792
 
 
d_1 =
 
 -3.0843959034513356629923904457422
  1.6374464759520504654890738027153
 

n =

        4963


omega_bar =

    0.4805


lambda_bar =

    0.9126


d_bar =

    1.1803

Example 3

Parameters as in Table 3 except for nu AND delta

nu = 15;          % capital to output ratio - MODIFIED
alpha = 0.025;    % productivity growth rate
beta = 0.02;      % population growth rate
delta = 0.10;     % depreciation rate - MODIFIED (NOTE: this is 0.05 in the paper)
r = 0.03;         % real interest rate
eta_p = 0.35;     % adjustment speed for prices
markup = 1.6;     % mark-up value
gamma = 0.8;      % (1-gamma) is the degree of money ilusion


% Numerical Values of Interest
c_eq=1-nu*(alpha+beta+delta)
eta_1 = fun_consumption_inv(c_eq)

A_w = eta_p*markup; % Coefficients in the quadratic equation A_w omega^2 + B_w omega + C_w = 0
B_w = alpha+beta-eta_p*(1+eta_1*markup);
C_w = eta_1*(eta_p-alpha-beta)+r*(eta_1-c_eq);

x=[];
syms x
omega_1 = vpasolve(A_w*x^2+B_w*x+C_w ==0,x)
i_omega_1 = eta_p*(markup*omega_1-1)
g_nominal_1 = alpha+beta+i_omega_1
lambda_1 = fun_Phi_inv(alpha+(1-gamma)*eta_p*(markup*omega_1-1))
d_1 = (omega_1-eta_1)./r

eta_1_calc=omega_1-r.*d_1; % verification step
c_eq_calc=fun_consumption(eta_1_calc); %verification step

% Sample paths of the dual Keen model

omega0 = 0.75;
lambda0 = 0.7;
d0 = 0.5;
Y0 = 100;
T = 200;

[tK,yK] = ode15s(@(t,y) dual_keen(y), [0 T], convert([omega0, lambda0, d0]), options);
yK = retrieve(yK);
Y_output = Y0*yK(:,2)/lambda0.*exp((alpha+beta)*tK);

n=length(yK(:,3))

omega_bar = yK(n,1)
lambda_bar = yK(n,2)
d_bar = yK(n,3)

g_bar=(1-c_minus)/nu-delta

figure(5)
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('dual_keen_phase_divergent','-depsc')


figure(6)
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('dual_keen_time_divergent','-depsc')
%

c_eq =

   -1.1750


eta_1 =

   0.0956 - 0.3934i

 
omega_1 =
 
  0.1375268793268625106771715577569 - 0.49617911095986802963355036671393i
 0.50275218583166006178002563156168 + 0.10281657191707506091722022599555i
 
 
i_omega_1 =
 
   - 0.27298494757695699402078392765613 - 0.2778603021375260965947882053598i
 - 0.06845877593427036540318564632546 + 0.057577280273562034113643326557506i
 
 
g_nominal_1 =
 
   - 0.22798494757695699402078392765613 - 0.2778603021375260965947882053598i
 - 0.02345877593427036540318564632546 + 0.057577280273562034113643326557506i
 
 
lambda_1 =
 
  0.97408289957607732385023306587867 - 0.021491278192165304543281057486094i
 0.96531792041637275646895074345848 + 0.0038394716601090056851942825208078i
 
 
d_1 =
 
 1.3963557103732356393956844389596 - 3.4272190639025027575057498472755i
 13.570532593866487342824153565785 + 16.539303698662266927519936576374i
 

n =

   830


omega_bar =

    0.0683


lambda_bar =

   4.7264e-08


d_bar =

   5.0107e+28


g_bar =

   -0.0353

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 y_h=omega-r*d disposable income of households

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, y]
    % It also work with the first 2 components 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) = old(:,1)-r.*old(:,3);
    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) = (old(:,1)-new(:,3))./r;
    end
end


function f = dual_keen(y)
    f = zeros(3,1);
    log_omega = y(1);
    tan_lambda = y(2);
    y_h = y(3);
    lambda = atan(tan_lambda)/pi+0.5;
    omega = exp(log_omega);
    g_Y = (1-fun_consumption(y_h))./nu-delta;

    f(1) = fun_Phi(lambda)-alpha-(1-gamma)*eta_p*(markup*omega-1); %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_consumption(y_h)-y_h)+(omega-y_h)*(g_Y+eta_p*(markup*omega-1)); %d(y_h)/dt
end

end