Implicit function in a for loop

Hi all,

I think I need your help. I have a model of three difference equations for three variables. Two initial values are needed, then the system is supposed to work through a for loop. Unfortunately, one of the main equations is implicit, so it must be solved numerically in every round. Here is what I have now:

clear
syms alpha beta delta theta A % parameters
alpha=1 % parameter values
beta=0.95
delta=0.1
theta=0.5
A=2
k1=10 % inital value for capital
h1=0.5 % initial value for labor supply
k(1)=k1 % 1st element of capital vector
h(1)=h1 % 1st element of labor supply vector
for t=1:100
    c(t)=((1-theta)/alpha)*(A^theta)*(k(t)^theta)/(h(t)^theta)*(1-h(t)) % endogenous consumption
    k(t+1)=(1-delta)*k(t)+(A^(1-theta))*((h(t))^(1-theta))*(k(t)^theta)-c(t)
    eqn=beta*(1-delta)+beta*theta*A^(1-theta)*(h(t+1)^(1-theta))/(k(t+1)^(1-theta))==(k(t+1)^theta)/(h(t+1)^theta)*(1-h(t+1))/(1-h(t))*(h(t)^theta)/(k(t)^theta)
    h(t+1)=solve(eqn,h(t+1))
end

Running ends up in this error message: Index exceeds the number of array elements. Index must not exceed 1.

How to proceed, guys? Any idea is warmly welcome. And thanks a lot.

Answer 1

[Walter Roberson]

syms alpha beta delta theta A % parameters
syms h_tplus1
Q = @(v) sym(v);
alpha = Q(1); % parameter values
beta = Q(0.95);
delta = Q(0.1);
theta = Q(0.5);
A = Q(2);
k1 = Q(10); % inital value for capital
h1 = Q(0.5); % initial value for labor supply
k(1) = double(k1); % 1st element of capital vector
h(1) = double(h1); % 1st element of labor supply vector
for t=1:100
    c(t) = double(((1-theta)/alpha)*(A^theta)*(k(t)^theta)/(h(t)^theta)*(1-h(t))); % endogenous consumption
    k(t+1) = double((1-delta)*k(t)+(A^(1-theta))*((h(t))^(1-theta))*(k(t)^theta)-c(t));
    eqn = beta*(1-delta)+beta*theta*A^(1-theta)*(h_tplus1^(1-theta))/(k(t+1)^(1-theta))==(k(t+1)^theta)/(h_tplus1^theta)*(1-h_tplus1)/(1-h(t))*(h(t)^theta)/(k(t)^theta);
    thissol = double(vpasolve(eqn, h_tplus1, [-inf inf]));  %constrain to real
    h(t+1) = thissol;
end
plot(h)

[Peter Galbacs]

Oh my goodness, it works. Many-many thanks. :)

[Peter Galbacs]

So I needed another solver - is it the problem?

[Walter Roberson]

I switched to vpasolve() for efficiency. solve() tries to find exact "algebraic" solutions -- closed form formulas.

The problem was that your eqn involved trying to use h(t+1) before you had assigned a value to h(t+1)

[Peter Galbacs]

Thanks a lot, Walter. An elegant solution. I will study your code to fully absorb the technical details. You helped a lot!

[Mendes Pereira]

The answer above does not represent the dynamics of h(t) in the original model. See further details below.

Answer 2

[Mendes Pereira]

The numerical simulation above is the simulation of something dynamic, but it is not the simulation of any element of the model described in the initial problem. If I am not making some mistake, the three equations in the original model are as follows:

Using the same parameter values as above, this system has a steady state given by: . However, the figure above shows that h(t) converges to 1 after a large number of iterations, while its actual steady state value is 0.426471.

I have a similar problem, but I do not know how to run a for loop with an implicit equation as part of the model. Some help would be very much appreciated.

[Torsten]

If h(t) converges to 1 like in the code from above, the behaviour of the term (1-h(t+1))/(1-h(t)) is unpredictable.