GARCH/asymm_garch_sim.f90 at main · Beliavsky/GARCH · GitHub

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
module asymm_garch_sim_mod
use        kind_mod, only: dp
use        util_mod, only: assert_equal, default
implicit none
private
public :: simulate_shift_twist_garch
integer, parameter :: nparam_gjr_garch = 5
contains

subroutine simulate_shift_twist_garch(z, ret, mu, omega, alpha, &
   gamma, beta, c_shift, sigma, iprint)
! simulate an asymmetric GARCH(1,1) process with shift and twist parameters
!
! Input:
! z     -- standardized noise
! mu    -- mean return
! omega -- variance constant
! alpha -- symmetric weight on past squared return
! gamma -- weight on past squared return when return is negative
! beta  -- weight on previous variance
! c_shift -- shift of standardized noise in variance equation
! iprint
!
! Output:
! ret -- Simulated return series
! sigma (optional) -- conditional standard deviation
real(kind=dp), intent(in)  :: mu, omega, alpha, gamma, beta, c_shift
real(kind=dp), intent(in)  :: z(:) ! noise with 0 mean and unit variance
real(kind=dp), intent(out) :: ret(:)
real(kind=dp), intent(out), optional :: sigma(:)
integer, intent(in), optional :: iprint
real(kind=dp) :: uc_var, sigma2, noise, noise_shift, pi_dp, phi_c, cdf_c, q
integer :: i, iprint_, n
real(kind=dp) :: sigma_i
iprint_ = default(0, iprint)
n = size(ret)
if (present(sigma)) call assert_equal(size(sigma), n, &
   "in simulate_shift_twist_garch, size(sigma)")
pi_dp = acos(-1.0_dp)
phi_c = exp(-0.5_dp*c_shift**2) / sqrt(2.0_dp*pi_dp)
cdf_c = 0.5_dp * (1.0_dp + erf(c_shift / sqrt(2.0_dp)))
q = (1.0_dp + c_shift**2) * cdf_c + c_shift * phi_c
uc_var = omega / (1.0_dp - beta - (1.0_dp + c_shift**2) * alpha - gamma * q) ! unconditional variance
if (iprint_ > 0) print*,"uc_var:", uc_var
sigma2 = uc_var
do i = 1, n
   sigma_i = sqrt(sigma2)
   noise = sigma_i * z(i)
   noise_shift = sigma_i * (z(i) - c_shift)
   if (present(sigma)) sigma(i) = sigma_i
   ret(i) = mu + noise
   if (iprint_ > 1) print "(i8, *(f10.6))", i,sigma_i, z(i), noise, ret(i)
   sigma2 = (alpha + merge(gamma,0.0_dp,noise_shift < 0.0_dp)) * &
              noise_shift**2 + beta*sigma2 + omega
end do
end subroutine simulate_shift_twist_garch
end module asymm_garch_sim_mod