Fortran to Matlab Translation with a dynamic fertility model

1 次查看（过去 30 天）
Sarah Schutz
Sarah Schutz 2019-12-14
回答： Ben Barrowes 2019-12-20
I am trying to translate fortran code to Matlab to solve a fertility problem using a dynamic discrete choiec model similar to Wolpin (1984). Here is an except of the Fortran code I am trying to translate. I have a copy of my Matlab code attached here but I cannot for the life of me figure out what I am doing wrong.
PROGRAM fertility
implicit none
external ran1,gasdev
integer numper,numvar,i,s,k,st1,st2,idnum,jj,j
integer calcsim
double precision alpha(6),vareps,pn
double precision summpi,lvalsum,lval1
double precision inc1,sr,c1,c2
! store shocks for emax calculation at each state point
! in each decision period (shocks are used for numerical
! integration) (need 9*100*20 = 18000 shocks)
double precision alleps(18000),epsvec(100)
double precision exval(9,20),utnc(100),utc(100)
double precision maxu(100),emaxu(0:8,1:20),sum1
! store shocks for simulation of behavior of
! 1000 women in 20 time periods (need 1000*20 periods=20000)
double precision simshock(20000)
double precision sutc,sutnc,simeps(20)
integer simchoice(1000,20),nn(1000,20)
double precision compfert
! initialize the parameters
alpha(1) = 0.000005
alpha(2) = 2000.0
alpha(3) = 0.0
alpha(4) = 0.0
alpha(5) = 0.0
alpha(6) = 2000.0
pn = 3000.0
alleps = 0.0
idnum = 15343 ! random number seed
vareps = (5000.0)**2
utnc = 0.0
utc = 0.0
! Draw set of shocks for each time period
call gasdev(idnum,alleps,18000)
! Solve the model backwards
! Will evaluate at all possible points
! in the state space, assuming the max number
! of children women can have is eight
oloop: do i = 20,1,-1
! loop over all possible state points (numbers of children)
sloop: do s = 0,8
! get residuals for each time period
! and possible state vector
st1 = (i-1)*800+s*100+1
st2 = st1+100-1
sr = s*1.0 ! convert he number of current children
! into a real number
epsvec = sqrt(vareps)*alleps(st1:st2)
inc1 = alpha(5)+alpha(6)*i ! get income income
c1 = inc1-pn*sr ! consumption if do not
! have another child
c2 = inc1-pn*(sr+1.0) ! consumption if do have
! another child
! calculate the utility for the set of 100 shocks if
! the woman has (utc) or does not have a child (utnc)
utnc = c1-0.5*alpha(1)*c1*c1+(alpha(2)+epsvec)*sr
utnc = utnc - alpha(3)*sr*sr+alpha(4)*c1*sr
utc = c2-0.5*alpha(1)*c2*c2+(alpha(2)+epsvec)*(sr+1.0)
utc = utc - alpha(3)*(sr+1.0)*(sr+1.0)+alpha(4)*c2*(sr+1.0)
! if not the terminal period, need to add in the emax
! values. If terminal period, only get current period
! utility. In either case, need to store the emax values
if (i.lt.20) then
utnc = utnc + emaxu(s,(i+1))
if (s.le.7) then
utc = utc + emaxu((s+1),(i+1))
else
utc = utnc ! if already have 8 children
end if ! cannot have any more
maxu = max(utnc,utc) ! calculate which choice gives
! the maximum utility
end if
if (i.eq.20) then ! if in terminal period, do not
if (s.eq.8) then ! need to add in future emax values
maxu = utnc
else
maxu = max(utnc,utc)
end if
end if
! calculate and store the emax values associated with each
! state point - the emax is determined by the mean of
! the utilities, taken over the realized shocks (which
! is numerically integrating with respect to the density of
! the shocks)
sum1 = 0.0
kloop: do k = 1,100 ! loop over the shocks
sum1 = sum1+maxu(k)
! write(*,*) ,epsvec(k),utc(k),utnc(k)
! if (maxu(k).eq.utnc(k)) then
! write(*,*) 'no child'
! else
! write(*,*) 'child'
! end if
end do kloop
emaxu(s,i) = sum1/100.0 ! calculate emax value
end do sloop
end do oloop
write(*,*) 'solved model - now doing simulation'

请先登录，再回答此问题。

回答（2 个）

A.B.
A.B. 2019-12-16
You did not tell us what error you get, so how do you expect people to be able to help you?
Ben Barrowes
Ben Barrowes 2019-12-20
You could try f2matlab here on the file exchange. You could also contact me at my author page if you still need help.
But you have cahnged the code substantially. Here is an untested conversion from f2matlab:
function fertility(varargin)
clear global;
clear functions;
global GlobInArgs nargs
GlobInArgs={mfilename,varargin{:}};
nargs=nargin+1;
persistent alleps alpha_fv c1 c2 calcsim compfert emaxu epsvec exval firstCall i idnum inc1 j jj k lval1 lvalsum maxu nn numper numvar pn s simchoice simeps simshock sr st1 st2 sum1 summpi sutc sutnc utc utnc vareps
;
if isempty(firstCall),firstCall=1;end;
if firstCall;
alleps=zeros(1,18000);
alpha_fv=zeros(1,6);
c1=0;
c2=0;
calcsim=0;
compfert=0;
emaxu=zeros(8+1,20);
epsvec=zeros(1,100);
exval=zeros(9,20);
i=0;
idnum=0;
inc1=0;
j=0;
jj=0;
k=0;
lval1=0;
lvalsum=0;
maxu=zeros(1,100);
nn=zeros(1000,20);
numper=0;
numvar=0;
pn=0;
s=0;
simchoice=zeros(1000,20);
simeps=zeros(1,20);
simshock=zeros(1,20000);
sr=0;
st1=0;
st2=0;
sum1=0;
summpi=0;
sutc=0;
sutnc=0;
utc=zeros(1,100);
utnc=zeros(1,100);
vareps=0;
end
% store shocks for emax calculation at each state point
% in each decision period (shocks are used for numerical
% integration) (need 9*100*20 = 18000 shocks)
% store shocks for simulation of behavior of
% 1000 women in 20 time periods (need 1000*20 periods=20000)
% initialize the parameters
firstCall=0;
alpha_fv(1) = 0.000005;
alpha_fv(2) = 2000.0;
alpha_fv(3) = 0.0;
alpha_fv(4) = 0.0;
alpha_fv(5) = 0.0;
alpha_fv(6) = 2000.0;
pn = 3000.0;
alleps(:) = 0.0;
% random number seed
idnum = 15343;
vareps =(5000.0).^2;
utnc(:) = 0.0;
utc(:) = 0.0;
% Draw set of shocks for each time period
[idnum,alleps]=gasdev(idnum,alleps,18000);
% Solve the model backwards
% Will evaluate at all possible points
% in the state space, assuming the max number
% of children women can have is eight
%oloop
for i = 20: -1: 1
% loop over all possible state points (numbers of children)
%sloop
for s = 0: 8
% get residuals for each time period
% and possible state vector
st1 =fix((i-1).*800 + s.*100 + 1);
st2 = fix(st1 + 100 - 1);
% convert he number of current children
sr = s.*1.0;
% into a real number
epsvec = sqrt(vareps).*alleps(st1:st2);
% get income income
inc1 = alpha_fv(5) + alpha_fv(6).*i;
% consumption if do not
c1 = inc1 - pn.*sr;
% have another child
% consumption if do have
c2 = inc1 - pn.*(sr+1.0);
% another child
% calculate the utility for the set of 100 shocks if
% the woman has (utc) or does not have a child (utnc)
utnc = c1 - 0.5.*alpha_fv(1).*c1.*c1 +(alpha_fv(2)+epsvec).*sr;
utnc = utnc - alpha_fv(3).*sr.*sr + alpha_fv(4).*c1.*sr;
utc = c2 - 0.5.*alpha_fv(1).*c2.*c2 +(alpha_fv(2)+epsvec).*(sr+1.0);
utc = utc - alpha_fv(3).*(sr+1.0).*(sr+1.0) + alpha_fv(4).*c2.*(sr+1.0);
% if not the terminal period, need to add in the emax
% values. If terminal period, only get current period
% utility. In either case, need to store the emax values
if(i < 20)
utnc = utnc + emaxu(s+1,(i+1));
if(s <= 7)
utc = utc + emaxu((s+1)+1,(i+1));
else;
% if already have 8 children
utc = utnc;
% cannot have any more
end
% calculate which choice gives
maxu = max(utnc,utc);
% the maximum utility
end
% if in terminal period, do not
if(i == 20)
% need to add in future emax values
if(s == 8)
maxu = utnc;
else;
maxu = max(utnc,utc);
end
end
% calculate and store the emax values associated with each
% state point - the emax is determined by the mean of
% the utilities, taken over the realized shocks (which
% is numerically integrating with respect to the density of
% the shocks)
sum1 = 0.0;
% loop over the shocks
%kloop
for k = 1: 100
sum1 = sum1 + maxu(k);
% write(*,*) ,epsvec(k),utc(k),utnc(k)
% if (maxu(k)==utnc(k)) then
% write(*,*) 'no child'
% else
% write(*,*) 'child'
% endif
%kloop
end
k = fix(100+1);
% calculate emax value
emaxu(s+1,i) = sum1./100.0;
%sloop
end
s = fix(8+1);
%oloop
clear all
end
i = fix(1+ -1);
[writeErrFlag]=writeFmt(1,['%c'],'''solved model - now doing simulation''');

类别

在 Help Center 和 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