## Contents

function out = inventories()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code simulates the Inventory Model as specified in
% Grasselli and Nguyen-Huu (2018)
%
%
% authors: M. Grasselli and A. Nguyen-Huu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc
warning off


## Base Parameters

This is a combination of Grasselli and Nguyen-Huu (2015) and Giraud and Grasselli (2019)

phi0 = 0.0401;       %Phillips curve parameter
phi1 = 6.41e-05;     %Phillips curve parameter

kappa0 = -0.0065;    % constant parameter in investment function
kappa1 = -5;         % affine parameter in exponent of investment function
kappa2 = 20;         % coefficient in the exponent of investment function
eta_u = 0.2;        % adjustment speed for investment as function of utilization
u_b = 0.9;           % baseline utilization

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;


## Other definitions

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


## Functions of interest

inflation = @(x1,x2,x3) eta_p*(markup*x1-1)+eta_q*(x2-x3);    % inflation, equation (31)

phillips = @(x) phi1./(1-x).^2- phi0;                         % Phillips curve, equation (32)
phillips_prime = @(x) 2*phi1./(1-x).^3;                       % derivative of phillips curve
phillips_inv = @(x) 1 - (phi1./(x+phi0)).^(1/2);              % inverse of phillips curve

investment= @(x1,x2) kappa0+exp(kappa1+kappa2.*x1)+eta_u*(x2-u_b); % investment function, equation (25)

% consumption, equation (24) with generalised logistic function depending
% on income only
consumption =@(x) max(c_minus,A_c + ((K_c-A_c)./(C_c+Q_c*exp(-B_c*x)).^(1/nu_c)));
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


## Example 1

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.8;     % adjustment speed for prices with respect to cost
eta_q = 0.2;     % adjustment speed for prices with respect to demand
eta_e = 0.8;     % adjustment speed for expected sales
eta_d = 0.7;     % adjustment speed for desired inventory
markup = 1.2;    % mark-up value
gamma = 0.4;     % (1-gamma) is the degree of money illusion
f_d = 0.1;       % desired inventories as a fraction of expected sales

% Functions of interest

inflation = @(x1,x2,x3) eta_p*(markup*x1-1)+eta_q*(x2-x3);    % inflation, equation (31)

phillips = @(x) phi1./(1-x).^2- phi0;                         % Phillips curve, equation (32)
phillips_prime = @(x) 2*phi1./(1-x).^3;                       % derivative of phillips curve
phillips_inv = @(x) 1 - (phi1./(x+phi0)).^(1/2);              % inverse of phillips curve

investment= @(x1,x2) kappa0+exp(kappa1+kappa2.*x1)+eta_u*(x2-u_b); % investment function, equation (25)

% consumption, equation (24) with generalised logistic function depending
% on income only
consumption =@(x) max(c_minus,A_c + ((K_c-A_c)./(C_c+Q_c*exp(-B_c*x)).^(1/nu_c)));
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

% Numerical Values of Interest

ydbar = 1/(1+(alpha+beta)*f_d);
yebar = 1/(1+(alpha+beta)*f_d)

x=fsolve(@system2d, [0.9,0.1]);

omegabar=x(1)
dbar=x(2)

ibar=inflation(omegabar,ydbar,yebar);

ubar=(nu*(alpha+beta+delta)*(1+(alpha+beta)*f_d))/(1-consumption(omegabar)*(1+(alpha+beta)*f_d))

kbar=investment(yebar*(1-omegabar)-r*dbar,ubar);

cbar = consumption(omegabar);

lambdabar = phillips_inv(alpha+(1-gamma)*ibar)

% Sample paths of the Keen model with inventories

omega0 = 0.9;
lambda0 = 0.9;
d0 = 0.5;
ye0 = 0.9;
u0 = 0.8;
y0=[convert([omega0, lambda0]),d0,ye0,u0];

Y0 = 100;
T = 200;

[tK,yK] = ode15s(@(t,y) keen_inventory(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);

n=length(yK(:,3));

omega_eq = yK(n,1)
lambda_eq = yK(n,2)
d_eq = yK(n,3)
ye_eq = yK(n,4)
u_eq = yK(n,5)

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(:,4),tK,yK(:,5))
ylabel('\omega,\lambda,y_e,u')

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

title(['\omega_0 = ', num2str(omega0,txt_format), ...
', \lambda_0 = ', num2str(lambda0, txt_format), ...
', ye_0 = ', num2str(ye0, txt_format), ...
', u_0 = ', num2str(u0, txt_format)])
legend('d','\omega', '\lambda','y_e','u')

% print('dual_keen_time_convergent','-depsc')

yebar =

0.9955

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.

omegabar =

0.8522

dbar =

0.3503

ubar =

1.8399

lambdabar =

0.9710

omega_eq =

0.8521

lambda_eq =

0.9710

d_eq =

0.3503

ye_eq =

0.9955

u_eq =

1.8399  ## 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 = keen_inventory(y)
f = zeros(5,1);
log_omega = y(1);
tan_lambda = y(2);
d = y(3);
ye = y(4);
u = y(5);

lambda = atan(tan_lambda)/pi+0.5;
omega = exp(log_omega);
pi_e = ye*(1-omega)-r*d;
yd = consumption(omega)+investment(pi_e,u)/u;
g_Y = (f_d*(alpha+beta+eta_d)+1)*(ye*(alpha+beta)+eta_e*(yd-ye))+eta_d*(yd-1);

f(1) = phillips(lambda)-alpha-(1-gamma)*inflation(omega,yd,ye); %d(log_omega)/dt
f(2) = (1+tan_lambda^2)*pi*lambda*(g_Y-alpha-beta); %d(tan_lambda)/dt
f(3) = d*(r-g_Y-inflation(omega,yd,ye))+omega-consumption(omega); %d(d)/dt
f(4) = ye*(alpha+beta-g_Y)+eta_e*(yd-ye);  % dy_e/dt
f(5) = u*(g_Y-investment(pi_e,u)/nu+delta);  % du/dt
end

function  F = system2d(x)
F(1) = investment(yebar*(1-x(1))-r*x(2),(nu*(alpha+beta+delta)*(1+(alpha+beta)*f_d))/(1-consumption(x(1))*(1+(alpha+beta)*f_d))) - nu*(alpha+beta+delta);
F(2) = x(2) - (x(1)-consumption(x(1)))/(alpha+beta+inflation(x(1),ydbar,yebar)-r);
end

end