This repository has been archived by the owner on Jan 30, 2023. It is now read-only.
-
-
Notifications
You must be signed in to change notification settings - Fork 7
/
markov_multifractal_cython.pyx
73 lines (56 loc) · 2.12 KB
/
markov_multifractal_cython.pyx
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
"""
Markov Switching Multifractal model
Cython code
"""
from sage.misc.randstate cimport randstate, current_randstate
cdef extern from "math.h":
double sqrt(double)
from .time_series cimport TimeSeries
def simulations(Py_ssize_t n, Py_ssize_t k,
double m0, double sigma,
int kbar, gamma):
"""
Return k simulations of length n using the Markov switching
multifractal model.
INPUT:
n, k -- positive integers
m0, sigma -- floats
kbar -- integer
gamma -- list of floats
OUTPUT:
list of lists
EXAMPLES::
sage: set_random_seed(0)
sage: msm = finance.MarkovSwitchingMultifractal(8,1.4,1.0,0.95,3)
sage: import sage.finance.markov_multifractal_cython
sage: sage.finance.markov_multifractal_cython.simulations(5,2,1.278,0.262,8,msm.gamma())
[[0.0014, -0.0023, -0.0028, -0.0030, -0.0019], [0.0020, -0.0020, 0.0034, -0.0010, -0.0004]]
"""
cdef double m1 = 2 - m0
cdef Py_ssize_t i, j, a, c
cdef TimeSeries t, eps
cdef TimeSeries markov_state_vector = TimeSeries(kbar)
cdef TimeSeries gamma_vals = TimeSeries(gamma)
cdef randstate rstate = current_randstate()
sigma = sigma / 100 # model's sigma is a percent
# output list of simulations
S = []
for i from 0 <= i < k:
# Initialize the model
for j from 0 <= j < kbar:
# n & 1 means "is odd"
markov_state_vector._values[j] = m0 if (rstate.c_random() & 1) else m1
t = TimeSeries(n)
# Generate n normally distributed random numbers with mean 0
# and variance 1.
eps = TimeSeries(n)
eps.randomize('normal')
for a from 0 <= a < n:
# Compute next step in the simulation
t._values[a] = sigma * eps._values[a] * sqrt(markov_state_vector.prod())
# Now update the volatility state vector
for c from 0 <= c < kbar:
if rstate.c_rand_double() <= gamma_vals._values[c]:
markov_state_vector._values[c] = m0 if (rstate.c_random() & 1) else m1
S.append(t)
return S