How to create a for-loop with matrices and vectors?

Dear Readers,

I am going to re-write this entire question. I have three equations describing steady state inflation, output and interest. The interest equation is an auxiliary equation, but still needed. My first step was to calculate al the steady state values for inflation and output and plot their relationship. The code used is as follows:

beta = 0.99 ; v = 21 ; gamma = 350 ; eps = 1 ; g = 0.2 ; sigma = 1;
c0 = 0.5 ;
pi = linspace(0.8,1.2,50) ;
for   i = 1 : length(pi)
pi_actual = pi(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, c0) ;
c0 = c_sol(i);
end
plot(pi,c_sol)

The code worked, but it somehow looks messy (?). And it had some errors within the third last line.

Now that I have all the solutions for c*, the steady state interest is calculated as pi* / beta. The next step is to look at the stability of all these 50 solutions by evaluating the Jacobian, an check if the eigenvalues are within the unit circle. Then I have to plot pi* versus these eigenvalues.

All help much appreciated!

PS.

Couldn't upload pictures anymore for next 24hrs.

[KSSV]

eig = eig(A); This is wrong, don't do it.....never name a variable on the name of a function.

[Guillaume]

Similarly, don't use pi as a variable name. It's also not clear if eval is a variable name or a call to the eval function. Whichever it is, don't do it.

[Stephen23]

Ouch!:

A = eval(J);
eig = eig(A);

https://blogs.mathworks.com/videos/2010/03/08/top-10-matlab-code-practices-that-make-me-cry/

@Sean Bagcik: this is not twitter, so why write ugly twitter-style tags? Just use a comma and no hash symbol.

[Cam Salzberger]

Fixed those tags for you.

Answer 1

[Cam Salzberger]

Hello Sean,

You can see the comments below your question for critiques on the code you posted. I would suggest a full revamp, but how you do it depends on how you are calculating your Jacobian. If you really wanted to vectorize the code, you could by popping out a 50*N x M matrix, lining up all the Jacobians on top of each other. If your Jacobian is a scalar (one function, one variable), then that's probably the best way.

Or you could change your three parameters into a cell array, and use cellfun to evaluate the whole thing. But for non-scalar Jacobians, I think sticking with a single for loop will be simpler, if not necessarily the fastest.

I would personally do it something like this:

% Assuming my parameters are x, y, and z
% Don't use i, j, pi, or any other built-in constants as variable names
x = rand(1,50); % For example
y = rand(1,50);
z = rand(1,50);

% Evaluate Jacobian, calculate eigen values, and store in a preallocated matrix
% Assuming the Jacobian will have 3 eigen values, for example
eigVals = zeros(3,numel(x));
for k = 1:numel(x) % Don't use i or j for variables
    J = createJacobian(x(k),y(k),z(k)); % You'll have to do your calculations for the Jacobian here
    eigVals(:,k) = eig(J);
end

% Plot against x
figure
axes
hold on
plot(x.',eigVals.') % Transposing because plot goes by column

Obviously this code won't run as is. You'll need to define the calculation for the Jacobian first.

Hope this helps!

-Cam

[Sean]

% Assuming my parameters are x, y, and z
x = linspace(0.8,1.05,50);
y = c_sol;
z = x/beta;

% Evaluate Jacobian, calculate eigen values, and store in a preallocated matrix
% Assuming the Jacobian will have 2 eigenvalues
eigVals = zeros(2,numel(x));
for k = 1:numel(x) % Don't use i or j for variables
    J = createJacobian(x(k),y(k),z(k)); % You'll have to do your calculations for the Jacobian here
    eigVals(:,k) = eig(J);
end

% Plot against x
figure
axes
hold on
plot(x.',eigVals.') % Transposing because plot goes by column

The 'createJacobian' and also 'jacobian' are not recognised when evaluating. What could be wrong here?

[Cam Salzberger]

Please notice this line in my original answer:

"Obviously this code won't run as is. You'll need to define the calculation for the Jacobian first."

I still don't fully understand your original question. However, if the function that you are looking to evaluate the Jacobian for is defined in your anonymous function, and if the only "variable" in it is pi_actual, then the Jacobian would just be a scalar at any given "pi_actual" value. Then the eigenvalue would be just one per, and it would be the same as the Jacobian evaluated at that point. So I kind of doubt that's the case.

So you'll need your own expression for how to calculate the Jacobian. My "createJacobian" function was just a placeholder.

[Sean]

beta = 0.99 ; v = 21 ; gamma = 350 ; eps = 1 ; g = 0.2 ; sigma = 1; alpha = 0.7;

c_sol=zeros(1,50);
x = linspace(0.8,1.2,50);
for   i = 1 : numel(x)
pi_actual = x(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
plot(x,c_sol)

z = x/beta;

This code gives me all the steady state solutions, 50 in total. The next step is to determine their stability. In order to do so I need to first evaluate a two-by-two Jacobian with symbolics (state vector is 2D) and after that implement all the solutions to determine the eigenvalues for each solution. The last step is plotting each inflation value versus the eigenvalue. The key is to show which inflation targets are 'stable'.

I find it difficult to calculate the Jacobian given the equations for inflation and output. Also I do not know how to do the steps after that.

Could you help here? Thanks!