CODEWITHSUNDEEP
matlabrecursion
Alex
Alex

Reputation: 39

Non-real root for a recursive difference equation in MATLAB

I was helped by user12339314 who helped me to solve a recursive difference equation in MATLAB. I have tried to apply the method to a slightly more complicated recursive equation, and it works too. This time, I have x(t+1) on the left-hand-side, and on the right-hand-side I have x(t) and x(t-1), i.e., the difference equation of order 3. And the method works, but I am little puzzled why I am getting a complex number, whose real part is the solution. I know the dynamics of xt are oscillatory, and it converges to 0.0480162880552655, but after I run the code below, and click on x, I see it converges to 0.0480162880552655 + 1.55851090378119e-67i. I am not sure why this happens. I tried different initial guesses, but the problem does not go away. Below is the code. Thank you for your help.

tic

clear
clc
format long;

time = 0:1:100;

A = 1.00;
h = 0.90;
beta  = 0.40;
alpha = 0.36;
n = 0.10;
xstar = 0.3581201041248495;
kstartheir = (alpha*A*xstar/h)^(1/(1-alpha));

zeta1 = 1/(1+beta+beta^2);
a = (1-zeta1+zeta1*beta^2*h/(1+n))/(zeta1+(1+h+n)*alpha/(h*(1-alpha)));
b = (zeta1*beta^2*h/(1+n))/(zeta1+(1+h+n)*alpha/(h*(1-alpha)));

recEq = @(p,q,z) ((1-alpha)*A*p^alpha/(1+h+n))*(1-zeta1*(1+h*q/(A*alpha*p^alpha))+(zeta1*beta^2*h/(1+n))*((1+h*p/(A*alpha*z^alpha))/(h*p/(A*alpha*z^alpha))));
x = nan(1, 100);
x(1) = 12200;
for t = 1:100
    kold = x(t);
    k = recEq(x(t),kold,kold);
    while abs(k-kold) > 1e-8
        kold = k;
        k = recEq(x(t), kold,kold);
    end
    x(t+1) = k;
end

toc

Upvotes: 0

Views: 64

Answers (1)

Thales
Thales

Reputation: 1316

Essentially, what you got, is a nonlinear recurrence equation. In order to get the x(t+1), you must solve (numerically) the recurrence equation to obtain it.

The idea user12339314 gave was a simple recurrence to find the fixed point of your equation. It is a simple root-finding algorithm.

Because now your recurrence equation has 3 terms, namely, x(t+1), x(t) and x(t-1), you must use these 3 values to find the root x(t+1). Changing a little your code, in order to make it clearer the "rood-finding part" of your code, we can write:

tic;clear;clc;format long;

time = 0:1:100;

A = 1.00;
h = 0.90;
beta  = 0.40;
alpha = 0.36;
n = 0.10;
xstar = 0.3581201041248495;
kstartheir = (alpha*A*xstar/h)^(1/(1-alpha));

zeta1 = 1/(1+beta+beta^2);
a = (1-zeta1+zeta1*beta^2*h/(1+n))/(zeta1+(1+h+n)*alpha/(h*(1-alpha)));
b = (zeta1*beta^2*h/(1+n))/(zeta1+(1+h+n)*alpha/(h*(1-alpha)));

recEq = @(p,q,z) ((1-alpha)*A*p^alpha/(1+h+n))*(1-zeta1*(1+h*q/(A*alpha*p^alpha))+(zeta1*beta^2*h/(1+n))*((1+h*p/(A*alpha*z^alpha))/(h*p/(A*alpha*z^alpha))));
x = nan(1, 100);
x(1) = 12200;
x(2) = 12200; % notice we must provide now 2 guesses, because we need 2 initial conditions

for t = 2:100
    k = x(t);
    kold = x(t);
    kold2 = x(t-1);

    % root finding algorithm to obtain x(t+1)
    converged = false;
    while ~converged
        kk = k;
        k = recEq(kold, kk, kold2);

        if(abs(k-kk)<1e-8)
            converged = true;
        end
    end
    x(t+1) = k;
end

toc

Notice a few things:

  • we must provide 2 initial conditions. So we define x(1) = 12200; and x(2) = 12200;
  • Given the equation you provided, you named in your recursive equation recEq the three variables: p=x(t), q=x(t+1) and z=x(t-1). Since you are looking for the root x(t+1), this is the variable k that you must define as the unknown. p and z are known (from the previous iterations or from the initial conditions in the first iteration).

Running this code, no round-off errors are introduced. The complex numbers are likely to appear in the exponentiation, with the wrong assignment (in your code) p,q,z in k = recEq(x(t), kold,kold);. Notice the difference from my code.

A last comment is that in the equation, you can use any root-finding algorithm. This is very interesting because, knowing that, you can use any built-in function to find the root. The advantage is that MATLAB's built-in functions are much better than any simple code you may write. We can change the loop to write it like:

for t = 2:100
    fun = @(k) recEq(x(t),k,x(t-1))-k;
    x(t+1) = fzero(fun,x(t));
    % x(t+1) = fsolve(fun,x(t));
end

To use the fzero or the fsolve built-in functions. You will get the same results (and for this particular case, fzero will run much faster than fsolve). If you have any other equation (even more difficult to solve), just use the built-in function and you should be fine.

Upvotes: 1

Related Questions