Fortran to Matlab Translation with a dynamic fertility model

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'

Answer 1

[A.B.]

You did not tell us what error you get, so how do you expect people to be able to help you?

Answer 2

[Ben Barrowes]

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''');