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
采纳的回答
Walter Roberson
2018-12-24
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.
6 个评论
Walter Roberson
2018-12-24
Hard to say. Your fortran code is completely invalid. It looks to me more like R code.
Walter Roberson
2018-12-24
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
2018-12-25
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
2018-12-25
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
2018-12-25
Might as well use the version I posted.
heidi pham
2018-12-25
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Fortran with MATLAB 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
