Is there a way to test the range of the parameters sigma and tau where my model remains stable?

I am attempting to find the range of values for the parameters sigma and tau where the model remains stable. Is there a way that this can be done in Matlab? Here is the code that I am working with:

tspan   = linspace(0, 100, 100001);
y0      =  0.7;   % initial value of state variable x1
r0      =  0.7;   % initial value of state variable x2
w0      =  0.8;   % initial value of state variable x3
x0      = [y0; r0; w0];
[t, x]  = ode45(@odefcn, tspan, x0);
%% Plot results
figure
subplot(3,1,1);
plot(t, x(:,1)); grid on
xlabel('Time'), ylabel('Output');
subplot(3,1,2);
plot(t, x(:,2)); grid on
xlabel('Time'), ylabel('Policy Rate');
subplot(3,1,3);
plot(t, x(:,3)); grid on
xlabel('Time'), ylabel('Wage Share');
dx  = zeros(numel(t), 3);
term  = zeros(numel(t), 1);
for i = 1:numel(t)
    [dx(i,:), term(i,:)]  = odefcn(t', x(i,:)');
end
%% Observe the effect of discontinuous term
figure
subplot(211)
plot(t, term),    grid on, title({'Effect of discontinuous term, $\delta y r + \tau y r$'}, 'interpreter', 'latex', 'fontsize', 12)
ylim([-1, 5])
subplot(212)
plot(t, dx(:,2)), grid on, title('dx_2')
ylim([-2, 4])
%% System of three differential equations
function [dx, term] = odefcn(t, x)
% definitions
y       = x(1);
r       = x(2);
w       = x(3);

% parameters
alpha   = 1.0;
beta    = 1.0;
gamma   = 0.6;
delta   = 0.6;
mu      = 0.1;
lambda  = 0.6;
theta   = 0.4;
omega   = 0.4;
sigma   = 0.1;
tau     = 0.1;

% ODEs
dx(1,1) = alpha * y - beta * r * y;    
dx(3,1) = - theta * w + lambda * y * w - mu * w * w;

if dx(1,1) <= 0 && dx(3,1) <= 0
    term    = 0;
else
    term    = delta * y * r + tau * y * r;
end

% Asymmetrical Reaction Function
if dx(1,1) <= 0 && dx(3,1) <= 0
    dx(2,1) = - omega * r + gamma * w * r + delta * y * r;
else
    % dx(2,1) = - omega * r + gamma * w * r + sigma * w * r + delta * y * r + tau * y * r; 
    dx(2,1) = - omega * r + gamma * w * r + sigma * w * r + term; 
end
end

[John D'Errico]

And, you know that it is a simple range for those parameters, and not some moderately complicated domain? I'd be a little surpised if it was something simple of course. But surely you know it is just a simple range for those parameters?

Instead, I would suggest you need to learn what causes that model to not be stable. Do a little reading, perhaps. Now look at the equations of state. Can you identify critical points where the process must be on the edge of stability?

[Christopher]

@John D'Errico Thanks for the comment. Of course! I just thought that there was something more expedient instead of assessing critical points (which I have done in my paper) and just manually changing the parmaters to assess the stability. In addition, I thought it might be something useful for more complex models. FWIW- I am just trying to do a robustness test to demonstrate the sensitivtiy of the model to the paramaters which is likely to be requested by journal referees.

[Torsten]

@Christopher

"Stable" means that you get a finite steady-state solution ?

Answer 1

[William Rose]

@Christopher,

One way to measure stability in your simulation is to compute the coefficient of variation (ratio of standard deviation to mean) for the second half of the simulation. I chose to do it for the second half to give the system time to stabilize. I chose to do it with state variable y=x(:,1)=output.

I added an outer pair of for loops to your script, to try different values of sigma and tau (sigma=0.00:0.05:0.50; tau=0.00:0.05:0.45.* I made a surface plot of the coefficient of variation, versus sigma and tau. I also plot the three state variables for the case of minimum and maximum coefficient of variation. See results below.

The surface plot above does not show clear regions of stability versus instability. It is more like a continuum.

Left plot: coeff. of variation for the second half of the simulation is maximum (sigma=0.25, tau=0.00).

Right plot: coeff. of var. for the second half of the simulation is minimum (sigma=0.45, tau=0.45).

If you don't like "coefficient of variation in the second half" as a stability measure, then make your own measure. You could still use the looping over sigma and tau, shown in this script, with a different measure of stability.

See attached script, and comments in it, for details. The script took about 18 seconds to run on my old notebook PC.

Good luck with your research.

* I use unequal number of values for sigma and tau, because when I make a surface plot of a square array, I sometimes get the axes reversed, without knowing it, and there is no error message. If the array is not square, then this silent error cannot happen.