Ticket #3356: sage-3356-part4.patch

sage/finance/markov_multifractal.py

@@ -26 +26 @@
             kbar   -- positive integer
             m0     -- float with 0 <= m0 <= 2
             sigma  -- positive float
-            gamma_kbar -- float with 0 < gamma_kbar < 1
+            gamma_kbar -- float with 0 <= gamma_kbar < 1
             b      -- float > 1
 
         EXAMPLES:
@@ -42 +42 @@
         self.__b = float(b)
         assert self.__b > 1, "b must be bigger than 1"
         self.__gamma_kbar = float(gamma_kbar)
-        assert self.__gamma_kbar > 0 and self.__gamma_kbar < 1, \
+        assert self.__gamma_kbar >= 0 and self.__gamma_kbar < 1, \
                "gamma_kbar must be between 0 and 1"
         self.__kbar = int(kbar)
         assert self.__kbar > 0, "kbar must be positive"
@@ -193 +193 @@
             Markov switching multifractal model with m0 = 1.278, sigma = 0.262, b = 2.11, and gamma_10 = 0.644
         """
         return markov_multifractal_cython.simulations(n, k,
-                   self.__m0, self.__sigma, self.__b,
+                   self.__m0, self.__sigma, 
                    self.__kbar, self.gamma())
         
 

sage/finance/markov_multifractal_cython.pyx

@@ -12 +12 @@
 from time_series cimport TimeSeries
 
 def simulations(Py_ssize_t n, Py_ssize_t k,
-               double m0, double sigma, double b,
+               double m0, double sigma, 
                int kbar, gamma):
     cdef double m1 = 2 - m0
     cdef Py_ssize_t i, j, a, c
@@ -21 +21 @@
     cdef TimeSeries gamma_vals = TimeSeries(gamma)
     cdef randstate rstate = current_randstate()
 
-    sigma = sigma / 100.0  # model's sigma is a percent
+    sigma = sigma / 100  # model's sigma is a percent
 
     # output list of simulations
     S = []
@@ -41 +41 @@
         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]:

sage/finance/time_series.pyx

@@ -1065 +1065 @@
         """
         return sqrt(self.variance(bias=bias))
 
+    def range_statistic(self):
+        sums = self.sums()
+        return (sums.max() - sums.min())/self.standard_deviation()
+
+    def hurst_exponent(self):
+        cdef double r = self.range_statistic()
+        if r == 0:
+            return 0
+        return log(r)/log(self._length)
+
     def min(self, bint index=False):
         """
         Return the smallest value in this time series. If this series
@@ -1243 +1253 @@
         for i from 0 <= i < bins:
             cnts[i] = 0
 
+        cdef Py_ssize_t j
         for i from 0 <= i < self._length:
-            cnts[<Py_ssize_t>floor((self._values[i] - mn)/step)] += 1
+            j = int((self._values[i] - mn)/step)
+            if j >= bins:
+                j = bins-1
+            cnts[j] += 1
 
         b = 1.0/(self._length * step)
         if normalize:
@@ -1276 +1290 @@
         s = 0
         for i, (x0,x1) in enumerate(intervals):
             s += polygon([(x0,0), (x0,counts[i]), (x1,counts[i]), (x1,0)], **kwds)
+        if len(intervals) > 0:
+            s.xmin(intervals[0][0])
+            s.xmax(intervals[-1][1])
+            s.ymin(0)
+            s.ymax(max(counts))
         return s
 
     def numpy(self):
@@ -1312 +1331 @@
         is platform independent.
         """
         if distribution == 'uniform':
-            self._randomize_uniform(loc, scale)
+            self._randomize_uniform(loc, loc + scale)
         elif distribution == 'normal':
             self._randomize_normal(loc, scale)
         elif distribution == 'semicircle':
@@ -1346 +1365 @@
         cdef Py_ssize_t k
         for k from 0 <= k < self._length:
             while 1:
-                x1 = 2.0 * rstate.c_rand_double() - 1.0
-                x2 = 2.0 * rstate.c_rand_double() - 1.0
-                w = x1 * x1 + x2 * x2
-                if w < 1.0: break
-            w = sqrt( (-2.0 * log( w ) ) / w )
+                x1 = 2*rstate.c_rand_double() - 1
+                x2 = 2*rstate.c_rand_double() - 1
+                w = x1*x1 + x2*x2
+                if w < 1: break
+            w = sqrt( (-2*log(w))/w )
             y1 = x1 * w
             y2 = x2 * w
             self._values[k] = m + y1*s