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')


%% 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')
% 

%% 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')
% 



%% 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