Ticket #3356: sage-3356-part2.patch

sage/finance/markov_multifractal.py

@@ -164 +164 @@
         self.__gamma = gamma
         return gamma
 
-    def simulation(self, T):
+    def simulation(self, n):
         """
-        Retun a time series of the T values taken by simulating the
-        running of this Markov switching multifractal model for T time
-        steps.
+        Same as self.simulations, but run only 1 time, and returns a time series
+        instead of a list of time series.
 
         INPUT:
-            T -- a positive integer
-        OUTPUT:
-            list -- a list of floats, the returns of the
-                    log-prices of a financial asset or
-                    exchange rate.
+            n -- a positive integer
 
         EXAMPLES:
-            sage: cad_usd = finance.MarkovSwitchingMultifractal(10,1.278,0.262,0.644,2.11); cad_usd
-            Markov switching multifractal model with m0 = 1.278, sigma = 0.262, b = 2.11, and gamma_10 = 0.644
-            sage: v = cad_usd.simulation(100)
-            sage: v    # random -- using seed doesn't work; planned rewrite of this function will work
-            [0.0011, -0.0032, 0.0006, 0.0007, 0.0034 ... -0.0023, 0.0008, 0.0015, -0.0003, 0.0027]
-            sage: v.sums()  # random
-            [0.0011, -0.0021, -0.0015, -0.0008, 0.0026 ... -0.0383, -0.0376, -0.0360, -0.0363, -0.0336]
-            sage: v.sums().exp().plot()
         """
-        # Two values of the distribution M.
-        m0 = self.m0()
-        vals = [m0, 2 - m0] 
+        return self.simulations(n, 1)[0]
 
-        # Initalize the Markov volatility state vector
-        from sage.rings.all import RDF
-        kbar = self.kbar()
-        m = (RDF**kbar)([random.choice(vals) for _ in xrange(kbar)])
-
-        sigma = self.sigma()/100.0
-
-        # r is the time series of returns
-        r = time_series.TimeSeries(T)
-
-        # Generate T Gaussian random numbers with mean 0
-        # and variance 1.
-        import scipy.stats
-        eps = scipy.stats.norm().rvs(T)
-
-        # Generate uniform distribution around 0 between -1 and 1
-        uniform = scipy.stats.uniform().rvs(kbar*T)
-
-        # The gamma_k
-        gamma = self.gamma()
-
-        for t in range(T):
-            r[t] = sigma * eps[t] * m.prod().sqrt()
-            # Now update the volatility state vector
-            j = t*kbar
-            for k in range(kbar):
-                if uniform[j+k] <= gamma[k]:
-                    # Draw from the distribution
-                    m[k] = random.choice(vals)
-                    
-        return r
-            
     def simulations(self, n, k=1):
         """
         Return k simulations of length n using this Markov switching
@@ -230 +183 @@
 
         INPUT:
             n -- positive integer; number of steps
-            k -- positive integer; number of simulations.
+            k -- positive integer (default: 1); number of simulations.
             
         OUTPUT:
             list -- a list of TimeSeries objects.
@@ -241 +194 @@
         """
         return markov_multifractal_cython.simulation(n, k,
                    self.__m0, self.__sigma, self.__b,
-                   self.__gamma_kbar, self.__kbar)
+                   self.__kbar, self.gamma())
         
 
     

sage/finance/markov_multifractal_cython.pyx

@@ -4 +4 @@
 Cython code
 """
 
-cdef extern from "stdlib.h":
-    long random()
+from sage.misc.randstate cimport randstate, current_randstate
+
+cdef extern from "math.h":
+    double sqrt(double)
 
 from time_series cimport TimeSeries
 
 def simulation(Py_ssize_t n, Py_ssize_t k,
                double m0, double sigma, double b,
-               double gamma_kbar, int kbar):
+               int kbar, gamma):
     cdef double m1 = 2 - m0
-    cdef Py_ssize_t i, j
-    cdef TimeSeries t
+    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.0  # model's sigma is a percent
 
@@ -25 +29 @@
     for i from 0 <= i < k:
         # Initalize the model
         for j from 0 <= j < kbar:
-            markov_state_vector._values[j] = m0 if random()%2 else m1
+            markov_state_vector._values[j] = m0 if (rstate.c_random() & 1) else m1    # n & 1 means "is odd"
         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
+            j = a * kbar
+            for c from 0 <= c < kbar:
+                if rstate.c_rand_double() <= gamma_vals._values[c]:
+                    markov_state_vector._values[k] = m0 if (rstate.c_random() & 1) else m1 
+
+        S.append(t)
+        
+    return S
+        
         
         
     

sage/finance/time_series.pyx

@@ -35 +35 @@
 include "../ext/cdefs.pxi"
 include "../ext/stdsage.pxi"
 include "../ext/python_string.pxi"
+include "../ext/random.pxi"
 
 cdef extern from "math.h":
     double exp(double)
@@ -45 +46 @@
 cdef extern from "string.h":
     void* memcpy(void* dst, void* src, size_t len)
 
-cdef extern from "stdlib.h":
-    long random()
-    long RAND_MAX
-
-cdef inline float ranf():
-    """
-    Return random number uniformly distributed in 0..1.
-    """
-    return (<float> random())/RAND_MAX
-    
 from sage.rings.integer import Integer
 from sage.rings.real_double import RDF
 from sage.modules.real_double_vector cimport RealDoubleVectorSpaceElement
@@ -916 +907 @@
         for i from 0 <= i < self._length:
             s += self._values[i]
         return s
+
+    def prod(self):
+        """
+        Return the prod of all the entries of self.  If self has
+        length 0, returns 1. 
+
+        OUTPUT:
+            double
+
+        EXAMPLES:
+            sage: v = finance.TimeSeries([1,1,1,2,3]); v
+            [1.0000, 1.0000, 1.0000, 2.0000, 3.0000]
+            sage: v.prod()
+            6.0
+        """
+        cdef double s = 1
+        cdef Py_ssize_t i
+        for i from 0 <= i < self._length:
+            s *= self._values[i]
+        return s
+        
     
     def mean(self):
         """
@@ -1051 +1063 @@
         else:
             return s
 
+    def clip(self, min=None, max=None):
+        """
+        Return new time series obtained from self by removing all values
+        <= a certain maximum value, >= a certain minimum, or both.
+
+        INPUT:
+            min -- None or double
+            max -- None or double
+
+        OUTPUT:
+            time seriessx
+        """
+        cdef Py_ssize_t i, j
+        cdef TimeSeries t
+        cdef double mn, mx
+        cdef double x
+        if min is None and max is None:
+            return self.copy()
+        # This code is ugly but I think as fast as possible.
+        if min is None and max is not None:
+            # Return everything <= max
+            j = 0
+            mx = max
+            for i from 0 <= i < self._length:
+                if self._values[i] <= mx:
+                    j += 1
+                    
+            t = TimeSeries(j)
+            j = 0
+            for i from 0 <= i < self._length:
+                if self._values[i] <= mx:
+                    t._values[j] = self._values[i]
+                    j += 1
+        elif max is None and min is not None:
+            # Return everything >= min
+            j = 0
+            mn = min
+            for i from 0 <= i < self._length:
+                if self._values[i] >= mn:
+                    j += 1
+            t = TimeSeries(j)
+            j = 0
+            for i from 0 <= i < self._length:
+                if self._values[i] >= mn:
+                    t._values[j] = self._values[i]
+                    j += 1
+        else:
+            # Return everything between min and max
+            j = 0
+            mn = min; mx = max
+            for i from 0 <= i < self._length:
+                x = self._values[i]
+                if x >= mn and x <= mx:
+                    j += 1
+            t = TimeSeries(j)
+            j = 0
+            for i from 0 <= i < self._length:
+                x = self._values[i]
+                if x >= mn and x <= mx:
+                    t._values[j] = x
+                    j += 1
+        return t
+
     def histogram(self, Py_ssize_t bins=50, bint normalize=False):
         """
         Return the frequency histogram of the values in
@@ -1118 +1193 @@
                 
         INPUT:
             bins -- positive integer (default: 50)
-            normalize -- whether to normalize so the total area in the
-                         bars of the histogram is 1. 
+            normalize -- (default: True) whether to normalize so the total
+                         area in the bars of the histogram is 1. 
             **kwds -- passed to the bar_chart function
         OUTPUT:
             a histogram plot
@@ -1153 +1228 @@
 
     def randomize(self, distribution='uniform', loc=0, scale=1, **kwds):
         """
+        Randomize the entries in this time series.
+        
         INPUT:
-            distribution -- 'uniform'
-                            'normal'
-                            'semicircle'
+            distribution -- 'uniform':    from loc to loc + scale
+                            'normal':     mean loc and standard deviation scale
+                            'semicircle': with center at loc (scale is ignored)
+            loc   -- float (default: 0)
+            scale -- float (default: 1)
+
+        NOTE: All random numbers are generated using the high quality
+        GMP random number function gmp_urandomb_ui.  This respects the
+        Sage set_random_state command.  It's not quite as fast at the
+        C library random number generator, but is better quality, and
+        is platform independent.
         """
         if distribution == 'uniform':
             self._randomize_uniform(loc, scale)
@@ -1170 +1255 @@
     def _randomize_uniform(self, double left, double right):
         if left >= right:
             raise ValueError, "left must be less than right"
+
+        cdef randstate rstate = current_randstate()
         cdef Py_ssize_t k
         cdef double d = right - left
         for k from 0 <= k < self._length:
-            self._values[k] = ranf() * d - left
+            self._values[k] = rstate.c_rand_double() * d + left
 
     def _randomize_normal(self, double m, double s):
         """
@@ -1185 +1272 @@
         """
         # Ported from http://users.tkk.fi/~nbeijar/soft/terrain/source_o2/boxmuller.c
         # This the box muller algorithm.
+        cdef randstate rstate = current_randstate()
+
         cdef double x1, x2, w, y1, y2
         cdef Py_ssize_t k
         for k from 0 <= k < self._length:
             while 1:
-                x1 = 2.0 * ranf() - 1.0
-                x2 = 2.0 * ranf() - 1.0
+                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 )
@@ -1204 +1293 @@
     def _randomize_semicircle(self, double center):
         cdef Py_ssize_t k
         cdef double x, y, s, d = 2, left = center - 1, z
+        cdef randstate rstate = current_randstate()
         z = d*d
         s = 1.5707963267948966192  # pi/2
         for k from 0 <= k < self._length:
             while 1:
-                x = ranf() * d - 1
-                y = ranf() * s
+                x = rstate.c_rand_double() * d - 1
+                y = rstate.c_rand_double() * s
                 if y*y + x*x < 1:
                     break
             self._values[k] = x + center