I have a system of three equations: inflation, output and interest. All variables are highly interconnected within these equations, however the interest equation is auxiliary.
The first step was determining the steady state combinations of inflation and consumption (see first image) and plotting that relationship. The following code was supposed to do so, but it gives an error about 'preallocation' of c_sol.
Above code should give me 50 combinations of inflation and output; interest is then calculated as inflation divided by beta.
The next step is evaluating the Jacobian at each of these values, determining the eigenvalues and plotting the eigenvalues with inflation on the x-axis. For stability it is required that the eigenvalues are within the unit circle, or in the case of complex numbers, the modulus. This is the code I have, but it needs a lot more editing:
Help is extremely appreciated!
https://i.stack.imgur.com/mNeKF.jpg
Update: this is the code I am working with, but I am missing the part after 'for'. I need to calculate a numerical matrix for each of those 50 values. 'c2' and 'p2' are lagged variables and equal to 'c' and 'p'.
syms c p R c2 p2
outp = (c2*p2)/(((99*R)/100 - 99/100)*(p2/p)^((2*R)/(R - 1))*(c2/c)^(R/(2*(R - 1))) + 99/100);
infl = (((2*c2*p2)/(35*(((99*R)/100 - 99/100)*(p2/p)^((2*R)/(R - 1))*(c2/c)^(R/(2*(R - 1))) + 99/100)) + 2/175)/(((99*R)/100 - 99/100)*(p2/p)^((2*R)/(R - 1))*(c2/c)^(R/(2*(R - 1))) + 99/100) - (99*p2)/100 + (99*p2^2)/100 + (3*((c2*p2)/(((99*R)/100 - 99/100)*(p2/p)^((2*R)/(R - 1))*(c2/c)^(R/(2*(R - 1))) + 99/100) + 1/5)^(20/7))/35 + 1/4)^(1/2) + 1/2;
J = [jacobian(outp, c2), jacobian(outp, p2) ; jacobian(infl, c2), jacobian(infl, p2)];
%%%%%%%%
beta = 0.99 ; v = 21 ; gamma = 350 ; eps = 1 ; g = 0.2 ; sigma = 1; alpha = 0.7;
c_sol = zeros(1,50);
p = linspace(0.8,1.2,50);
for i = 1 : numel(p)
pi_actual = p(i);
fun = @(c) (1 - beta)*pi_actual*(pi_actual-1) - (v/(alpha*gamma))*(c+g).^((1+eps)/alpha) - ((1-v)/gamma)*(c+g)*c.^(-sigma);
c_sol(i) = fzero(fun,0.5);
end
c = c_sol; c2 = c_sol; p2 = p;
R = p/beta; % interest
%%%%%%%
eigVals = zeros(2,numel(p));
for k = 1:numel(p)
eigVals(:,k) = eig(B);
end