-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathps2.f95
122 lines (97 loc) · 4.33 KB
/
ps2.f95
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
! Last change: PND 6 Sep 2007 5:56 pm
module parameters
implicit none
REAL, PARAMETER :: b = 0.99, d = 0.025, a = 0.36
REAL, PARAMETER :: klb = 0.01, inc = 0.025, kub = 45.0
INTEGER, PARAMETER :: length_grid_k = (kub-klb)/inc+1
INTEGER, PARAMETER:: no_states = 2.
REAL, PARAMETER:: stay_good_prob = 0.977, stay_bad_prob = 0.926 ! for the transitions
REAL, PARAMETER :: toler = 1.e-4 ! Numerical tolerance
end module
! ============================================================================================
! ============================================================================================
PROGRAM HW2NonStochastic
REAL :: total, dist
REAL, DIMENSION(2) :: elapsed
call solution
call etime(elapsed, total)
PRINT*,'--------------------------------------------------'
PRINT*,'total time elapsed =',elapsed, total
PRINT*,'--------------------------------------------------'
END PROGRAM HW2NonStochastic
! ============================================================================================
subroutine solution
USE parameters
IMPLICIT NONE
REAL, dimension(:), allocatable :: Kgrid, g_k
REAL, dimension(:,:), allocatable:: value, value_new
INTEGER, DIMENSION(:,:),allocatable :: PFUNC
REAL :: Zgrid(no_states) = (/1.25, 0.2/)
REAL :: Pi(no_states, no_states)
REAL, dimension(:,:,:), allocatable :: vtmp
INTEGER:: iter, index_k, index_kp, index_z
REAL:: diffg, diffb, diff, k, kp, c, z
INTEGER:: i = 1
allocate(value(length_grid_k, no_states))
allocate(value_new(length_grid_k, no_states))
allocate(PFUNC(length_grid_k, no_states))
allocate(vtmp(length_grid_k, length_grid_k, no_states))
allocate(Kgrid(length_grid_k))
allocate(g_k, mold = Kgrid)
do while (i <= length_grid_k) !do loop for assigning capital grid K
Kgrid(i) = klb + (i-1)*inc
!write(*,*) i, Kgrid(i)
i = i+1
end do
Pi(1, :) = (/stay_good_prob, 1-stay_bad_prob/)
Pi(2, :) = (/1-stay_good_prob, stay_bad_prob/)
iter = 1
diff = 1000.d0
value(:,1)= 0.*Kgrid !Initial Value guess
value(:, 2) = 0.*Kgrid
do while (diff >= toler)
!$omp parallel do default(shared) private(index_k, index_z, k, vtmp, z, kp, c)
do index_k = 1, length_grid_k ! Capital grid
k = Kgrid(index_k)
vtmp(index_k, :,:) = -1.0e-16
do index_z = 1, no_states
z = Zgrid(index_z)
do index_kp = 1, length_grid_k
kp = Kgrid(index_kp)
c = k**a+(1.-d)*k-kp
if (c > 0.) then
vtmp(index_k, index_kp, index_z) = log(c)+b*(Pi(index_z, 1)*value(index_kp, 1)+Pi(index_z, 2)*value(index_kp, 2))
endif
enddo
value_new(index_k, index_z) = MAXVAL(vtmp(index_k, :,index_z))
PFUNC(index_k, index_z) = MAXLOC(vtmp(index_k, :, index_z), 1)
g_k(index_k) = Kgrid(MAXLOC(vtmp(index_k, :, index_z), 1))
enddo
enddo
!$omp end parallel do
diffg = maxval(abs(value_new(:,1)-value(:,1)))/ABS(value_new(length_grid_k, 1))
diffb = maxval(abs(value_new(:,2)-value(:,2)))/ABS(value_new(length_grid_k, 2))
diff = max(abs(diffg), abs(diffb))
value = value_new
iter = iter+1
enddo
print *, ' '
print *, 'Successfully converged with sup_norm ', diff, 'with number of iteration:', iter
!print *, g_k
!CALL vcDrawCurve(d, Kgrid, g_k, length_grid_k)
open (UNIT = 1, FILE='VFUNC',STATUS='replace')
do index_k = 1, length_grid_k
WRITE(UNIT = 1, FMT=*) value(index_k, :)
end do
close (UNIT = 1)
open (UNIT = 2, FILE='AGRID',STATUS='replace')
do index_k = 1, length_grid_k
WRITE(UNIT = 2, FMT=*) Kgrid(index_k)
end do
close (UNIT = 2)
open (UNIT = 3, FILE='PFUNC',STATUS='replace')
do index_k = 1, length_grid_k
WRITE(UNIT = 3, FMT=*) PFUNC(index_k,:)
end do
close (UNIT = 3)
end subroutine