! ---------------------------------------------------------------------- ! The routines for random number generation were obtained from ! the book NUMERICAL RECIPES SUBROUTINE ran1sub(idum,x) IMPLICIT NONE INTEGER, PARAMETER :: K4B=selected_int_kind(9) INTEGER(K4B), INTENT(INOUT) :: idum REAL(8) :: x INTEGER(K4B), PARAMETER :: IA=16807, IM=2147483647,IQ=127773,IR=2836 REAL(8), SAVE :: am INTEGER(K4B), SAVE :: ix=-1,iy=-1,k if (idum <= 0 .or. iy < 0) then am=nearest(1.0,-1.0)/IM iy=ior(ieor(888889999,abs(idum)),1) ix=ieor(777755555,abs(idum)) idum=abs(idum)+1 end if ix=ieor(ix,ishft(ix,13)) ix=ieor(ix,ishft(ix,-17)) ix=ieor(ix,ishft(ix,5)) k=iy/IQ iy=IA*(iy-k*IQ)-IR*k if (iy<0) iy=iy+IM x=am*ior(iand(IM,ieor(ix,iy)),1) END SUBROUTINE ran1sub SUBROUTINE ran1(idumtemp,vv,n) implicit none external ran1sub REAL(8) vv(n),x INTEGER idumtemp,kk,n vvloop: do kk = 1,n call ran1sub(idumtemp,x) vv(kk) = x end do vvloop END SUBROUTINE ran1 !------------------------ SUBROUTINE GASDEV(idum,vv,n) ! This function returns a normally distributed deviate with ! zero mean and unit variance, using ran1(IDUM) as the ! source of uniform deviates integer iset,n,k,gotoc real(8) v1,v2,r,fac,gset,gasdev1 real(8) temp1(2),vv(n) external ran1 write(*,*) 'want',n,'random numbers' vvloop: do k = 1,n iset = 0 v1 = 0. v2 = 0. r = 0. fac = 0. gset = 0. gasdev1 = 0. ix1=0. ix2=0. gotoc = 0 if (iset.eq.0) then 1 call ran1(idum,temp1,2) v1 = 2.*temp1(1)-1. v2 = 2.*temp1(2)-1. ! write(*,*) 'v1',v1 ! write(*,*) 'v2',v2 ! write(*,*) 'idum',idum r = v1**2.+v2**2. if (r.ge.1.) then gotoc = gotoc + 1 if (gotoc.gt.n) then write(*,*) 'error in gasdev' end if go to 1 end if fac = dsqrt(-2.*dlog(r)/r) gset = v1*fac gasdev1 = v2*fac iset = 1 else gasdev1 = gset iset = 0 end if vv(k)=gasdev1 ! write(*,*) 'random normal number',k,vv(k) end do vvloop return end SUBROUTINE GASDEV !---------------------------------------------------------------------------- ! This program implements a dynamic fertility choice model. ! There are twenty periods (each representing two years) and ! women choose in each period whether to have a child or not. ! They can have up to eight children. 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' !******************************************************************************* ! Simulation part of the program ! ! The above calculation yields the emax values for all the time periods. ! Now we can use them to simulate the behavior of 1000 women under ! the model: calcsim = 1 if (calcsim.eq.1) then call gasdev(idnum,simshock,20000) ! draw shocks simchoice = 0.0 simloop: do i =1,1000 st1 = (i-1)*20+1 st2 = st1+20-1 simeps=sqrt(vareps)*simshock(st1:st2) yrloop: do j = 1,20 ! calculate number of children based on ! previous choices if (j.eq.1) then sr = 0.0 else sr = sum(simchoice(i,1:(j-1))) end if inc1 = alpha(5)+alpha(6)*j ! income c1 = inc1-pn*sr ! consumption c2 = inc1-pn*(sr+1.0) ! calculate the utility for the current period shock if ! the woman has (sutc) or does not have a child (sutnc) sutnc = c1-0.5*alpha(1)*c1*c1+(alpha(2)+simeps(j))*sr sutnc = sutnc - alpha(3)*sr*sr+alpha(4)*c1*sr if (j.lt.20) then sutnc = sutnc + emaxu(sr,j+1) end if simchoice(i,j) = 0 if (sr.le.7.0) then sutc = c2-0.5*alpha(1)*c2*c2 sutc = sutc + (alpha(2)+simeps(j))*(sr+1.0) sutc = sutc - alpha(3)*(sr+1.0)*(sr+1.0) sutc = sutc + alpha(4)*c2*(sr+1.0) if (j.lt.20) then ! add in emax if not in terminal period sutc = sutc + emaxu((sr+1),j+1) end if if (sutc.gt.sutnc) then simchoice(i,j)= 1 end if end if end do yrloop end do simloop end if ! calcsim ! write out the simulated choices for 10 of the women rloop: do i=1,10 write(*,*) 'person',i write(*,*) (simchoice(i,j),j=1,20) write(*,*) 'total number of children' write(*,*) sum(simchoice(i,:)) end do rloop ! calculate summary statistics - average number of children and average number ! by age of the woman rrloop: do i = 1,1000 zzloop: do j = 1,20 nn(i,j) = sum(simchoice(i,1:j)) end do zzloop end do rrloop write(*,*) 'Avg fertility by model period' ggloop: do i = 1,20 compfert = sum(nn(:,i))/1000.00 write(*,15) i, compfert end do ggloop 15 format(i5,f8.2) end program fertility