Convert fortran code to Matlab code C - Get result as complex numbers

Hello,

I tried to convert a fortran code to matlab code. I don't know why I get very weird resutl after running my Malab code (the result comprise of complex numbers).

Here is the fortran code http://www.jonathanheathcote.com/Macro2/FYMHW1

If someone could see what did I do wrong in converting, please kindly help me. Thanks so much.

%  set parameter values
%
theta  = 0.4;               % production parameter
beta   = 0.95;              % subjective discount factor 
delta  = 0.1;               % capital depreciation rate
A=6.0476;

% Horizon
T=38; %which is "capt" in fortran code

% Sequence of capital and consumption
kseq=zeros(T+1,1);
cseq=zeros(T+1,1);

%Initial capital
kseq(1)=10;

% Guess the range of initial consumption
cguessl=1;
cguessh=10;

%Iteration parameter
iter=0;
tol=0.01;

while gapk > 0.01 && iter < 1000
    cguess=0.5*(cguessl+cguessh);
    cseq(1)=cguess;
   for i = 2:T+1
        %Law of motion of capital
        kseq(i) = A*(kseq(i-1))^theta+(1-delta)*kseq(i-1)-cseq(i-1);
        
        %check whether capital goes negative, if so reduce guess for
        %consumption c(1)
        if kseq(i)<0
            cguessh=cguess;            
        else
        cseq(i)=beta*(1+theta*A*(kseq(i))^(theta-1)-delta)*cseq(i-1);
        end
   end
   %Check terminal capital and adjust bounds of cguess effectively
   if kseq(T+1)>0
       cguessl=cguess;
   else
       cguessh=cguess;
   end
          gapk = abs(kseq(T+1));
          iter=iter+1;
end
   

Answer 1

[Walter Roberson]

You have

        kseq(i) = (kseq(i-1))^theta+(1-delta)*kseq(i-1)-cseq(i-1);
        
        %check whether capital goes negative, if so reduce guess for
        %consumption c(1)
        if kseq(i)<0
            cguessh=cguess;            
        else
        cseq(i)=beta*(1+theta*(kseq(i))^(theta-1)-delta)*cseq(i-1);
        end

In the case kseq(i) < 0, you set cguessh but you do not change kseq(i) . Then on the next iteration of the for i, kseq(i-1) will refer to the negative value, and you raise the negative value to theta. MATLAB defines the element-wise ^ operator as A.^B = exp(log(A)*B) and so when A < 0 the log(A) is complex and the whole result ends up complex unless B is an integer, even if there is a negative real number C such that C^(1/B) = A.

[Walter Roberson]

Hard to say. Your fortran code is completely invalid. It looks to me more like R code.

[Walter Roberson]

   if kseq(T+1)>0; flag=2;

You are doing that test even if you broke the for i early, in which case you did not update kseq(T+1)

Note that your if/elseif sequence assumes that the final else can only occur if kseq(T+1) equals 0, but that branch can be reached if kseq(T+1) becomes nan.

[Walter Roberson]

kseq(T+1)>0 can not be true if kseq(T+1) is nan.

You can test for nan by using isnan(kseq(T+1))

but your bigger problem is that you are testing kseq(T+1) even when you did not change kseq(T+1) because you exited the loop early.

[Walter Roberson]

while gapk > 0.01 && iter < 1000
    cguess=0.5*(cguessl+cguessh);
    cseq(1)=cguess;
    flag = 0;
    for i = 2:T+1
        %Law of motion of capital
        kseq(i) = A*(kseq(i-1))^theta+(1-delta)*kseq(i-1)-cseq(i-1);
        %check whether capital goes negative, if so reduce guess for
        %consumption c(1)
        if kseq(i)<0; flag=1; break; end
        cseq(i)=beta*(1+theta*A*(kseq(i))^(theta-1)-delta)*cseq(i-1);
    end
    %Check terminal capital and adjust bounds of cguess effectively
    if flag == 0 && kseq(T+1)>0
        cguessl = cguess;
    else
        cguessh = cguess;
    end

    gapk = abs(kseq(T+1));
    iter=iter+1;
end

[Walter Roberson]

Might as well use the version I posted.

[heidi pham]

thank you!