Ticket #3356: sage-3356-part2.patch

File sage-3356-part2.patch, 13.0 KB (added by was, 14 years ago)
  • sage/finance/markov_multifractal.py

    # HG changeset patch
    # User William Stein <wstein@gmail.com>
    # Date 1212514292 25200
    # Node ID 1f137599f3e3141f261529ab6930ae9880b2443d
    # Parent  b76f5b58a6addd3c0a302b79d7103d90aed4a435
    trac #3356 finance -- more work on quality random number generation; vastly optimized multifractal model
    
    diff -r b76f5b58a6ad -r 1f137599f3e3 sage/finance/markov_multifractal.py
    a b class MarkovSwitchingMultifractal: 
    164164        self.__gamma = gamma
    165165        return gamma
    166166
    167     def simulation(self, T):
     167    def simulation(self, n):
    168168        """
    169         Retun a time series of the T values taken by simulating the
    170         running of this Markov switching multifractal model for T time
    171         steps.
     169        Same as self.simulations, but run only 1 time, and returns a time series
     170        instead of a list of time series.
    172171
    173172        INPUT:
    174             T -- a positive integer
    175         OUTPUT:
    176             list -- a list of floats, the returns of the
    177                     log-prices of a financial asset or
    178                     exchange rate.
     173            n -- a positive integer
    179174
    180175        EXAMPLES:
    181             sage: cad_usd = finance.MarkovSwitchingMultifractal(10,1.278,0.262,0.644,2.11); cad_usd
    182             Markov switching multifractal model with m0 = 1.278, sigma = 0.262, b = 2.11, and gamma_10 = 0.644
    183             sage: v = cad_usd.simulation(100)
    184             sage: v    # random -- using seed doesn't work; planned rewrite of this function will work
    185             [0.0011, -0.0032, 0.0006, 0.0007, 0.0034 ... -0.0023, 0.0008, 0.0015, -0.0003, 0.0027]
    186             sage: v.sums()  # random
    187             [0.0011, -0.0021, -0.0015, -0.0008, 0.0026 ... -0.0383, -0.0376, -0.0360, -0.0363, -0.0336]
    188             sage: v.sums().exp().plot()
    189176        """
    190         # Two values of the distribution M.
    191         m0 = self.m0()
    192         vals = [m0, 2 - m0]
     177        return self.simulations(n, 1)[0]
    193178
    194         # Initalize the Markov volatility state vector
    195         from sage.rings.all import RDF
    196         kbar = self.kbar()
    197         m = (RDF**kbar)([random.choice(vals) for _ in xrange(kbar)])
    198 
    199         sigma = self.sigma()/100.0
    200 
    201         # r is the time series of returns
    202         r = time_series.TimeSeries(T)
    203 
    204         # Generate T Gaussian random numbers with mean 0
    205         # and variance 1.
    206         import scipy.stats
    207         eps = scipy.stats.norm().rvs(T)
    208 
    209         # Generate uniform distribution around 0 between -1 and 1
    210         uniform = scipy.stats.uniform().rvs(kbar*T)
    211 
    212         # The gamma_k
    213         gamma = self.gamma()
    214 
    215         for t in range(T):
    216             r[t] = sigma * eps[t] * m.prod().sqrt()
    217             # Now update the volatility state vector
    218             j = t*kbar
    219             for k in range(kbar):
    220                 if uniform[j+k] <= gamma[k]:
    221                     # Draw from the distribution
    222                     m[k] = random.choice(vals)
    223                    
    224         return r
    225            
    226179    def simulations(self, n, k=1):
    227180        """
    228181        Return k simulations of length n using this Markov switching
    class MarkovSwitchingMultifractal: 
    230183
    231184        INPUT:
    232185            n -- positive integer; number of steps
    233             k -- positive integer; number of simulations.
     186            k -- positive integer (default: 1); number of simulations.
    234187           
    235188        OUTPUT:
    236189            list -- a list of TimeSeries objects.
    class MarkovSwitchingMultifractal: 
    241194        """
    242195        return markov_multifractal_cython.simulation(n, k,
    243196                   self.__m0, self.__sigma, self.__b,
    244                    self.__gamma_kbar, self.__kbar)
     197                   self.__kbar, self.gamma())
    245198       
    246199
    247200   
  • sage/finance/markov_multifractal_cython.pyx

    diff -r b76f5b58a6ad -r 1f137599f3e3 sage/finance/markov_multifractal_cython.pyx
    a b Cython code 
    44Cython code
    55"""
    66
    7 cdef extern from "stdlib.h":
    8     long random()
     7from sage.misc.randstate cimport randstate, current_randstate
     8
     9cdef extern from "math.h":
     10    double sqrt(double)
    911
    1012from time_series cimport TimeSeries
    1113
    1214def simulation(Py_ssize_t n, Py_ssize_t k,
    1315               double m0, double sigma, double b,
    14                double gamma_kbar, int kbar):
     16               int kbar, gamma):
    1517    cdef double m1 = 2 - m0
    16     cdef Py_ssize_t i, j
    17     cdef TimeSeries t
     18    cdef Py_ssize_t i, j, a, c
     19    cdef TimeSeries t, eps
    1820    cdef TimeSeries markov_state_vector = TimeSeries(kbar)
     21    cdef TimeSeries gamma_vals = TimeSeries(gamma)
     22    cdef randstate rstate = current_randstate()
    1923
    2024    sigma = sigma / 100.0  # model's sigma is a percent
    2125
    def simulation(Py_ssize_t n, Py_ssize_t  
    2529    for i from 0 <= i < k:
    2630        # Initalize the model
    2731        for j from 0 <= j < kbar:
    28             markov_state_vector._values[j] = m0 if random()%2 else m1
     32            markov_state_vector._values[j] = m0 if (rstate.c_random() & 1) else m1    # n & 1 means "is odd"
    2933        t = TimeSeries(n)
     34       
     35        # Generate n normally distributed random numbers with mean 0
     36        # and variance 1.
     37        eps = TimeSeries(n)
     38        eps.randomize('normal')
     39
     40        for a from 0 <= a < n:
     41            # Compute next step in the simulation
     42            t._values[a] = sigma * eps._values[a] * sqrt(markov_state_vector.prod())
     43
     44            # Now update the volatility state vector
     45            j = a * kbar
     46            for c from 0 <= c < kbar:
     47                if rstate.c_rand_double() <= gamma_vals._values[c]:
     48                    markov_state_vector._values[k] = m0 if (rstate.c_random() & 1) else m1
     49
     50        S.append(t)
     51       
     52    return S
     53       
    3054       
    3155       
    3256   
  • sage/finance/time_series.pyx

    diff -r b76f5b58a6ad -r 1f137599f3e3 sage/finance/time_series.pyx
    a b include "../ext/cdefs.pxi" 
    3535include "../ext/cdefs.pxi"
    3636include "../ext/stdsage.pxi"
    3737include "../ext/python_string.pxi"
     38include "../ext/random.pxi"
    3839
    3940cdef extern from "math.h":
    4041    double exp(double)
    cdef extern from "string.h": 
    4546cdef extern from "string.h":
    4647    void* memcpy(void* dst, void* src, size_t len)
    4748
    48 cdef extern from "stdlib.h":
    49     long random()
    50     long RAND_MAX
    51 
    52 cdef inline float ranf():
    53     """
    54     Return random number uniformly distributed in 0..1.
    55     """
    56     return (<float> random())/RAND_MAX
    57    
    5849from sage.rings.integer import Integer
    5950from sage.rings.real_double import RDF
    6051from sage.modules.real_double_vector cimport RealDoubleVectorSpaceElement
    cdef class TimeSeries: 
    916907        for i from 0 <= i < self._length:
    917908            s += self._values[i]
    918909        return s
     910
     911    def prod(self):
     912        """
     913        Return the prod of all the entries of self.  If self has
     914        length 0, returns 1.
     915
     916        OUTPUT:
     917            double
     918
     919        EXAMPLES:
     920            sage: v = finance.TimeSeries([1,1,1,2,3]); v
     921            [1.0000, 1.0000, 1.0000, 2.0000, 3.0000]
     922            sage: v.prod()
     923            6.0
     924        """
     925        cdef double s = 1
     926        cdef Py_ssize_t i
     927        for i from 0 <= i < self._length:
     928            s *= self._values[i]
     929        return s
     930       
    919931   
    920932    def mean(self):
    921933        """
    cdef class TimeSeries: 
    10511063        else:
    10521064            return s
    10531065
     1066    def clip(self, min=None, max=None):
     1067        """
     1068        Return new time series obtained from self by removing all values
     1069        <= a certain maximum value, >= a certain minimum, or both.
     1070
     1071        INPUT:
     1072            min -- None or double
     1073            max -- None or double
     1074
     1075        OUTPUT:
     1076            time seriessx
     1077        """
     1078        cdef Py_ssize_t i, j
     1079        cdef TimeSeries t
     1080        cdef double mn, mx
     1081        cdef double x
     1082        if min is None and max is None:
     1083            return self.copy()
     1084        # This code is ugly but I think as fast as possible.
     1085        if min is None and max is not None:
     1086            # Return everything <= max
     1087            j = 0
     1088            mx = max
     1089            for i from 0 <= i < self._length:
     1090                if self._values[i] <= mx:
     1091                    j += 1
     1092                   
     1093            t = TimeSeries(j)
     1094            j = 0
     1095            for i from 0 <= i < self._length:
     1096                if self._values[i] <= mx:
     1097                    t._values[j] = self._values[i]
     1098                    j += 1
     1099        elif max is None and min is not None:
     1100            # Return everything >= min
     1101            j = 0
     1102            mn = min
     1103            for i from 0 <= i < self._length:
     1104                if self._values[i] >= mn:
     1105                    j += 1
     1106            t = TimeSeries(j)
     1107            j = 0
     1108            for i from 0 <= i < self._length:
     1109                if self._values[i] >= mn:
     1110                    t._values[j] = self._values[i]
     1111                    j += 1
     1112        else:
     1113            # Return everything between min and max
     1114            j = 0
     1115            mn = min; mx = max
     1116            for i from 0 <= i < self._length:
     1117                x = self._values[i]
     1118                if x >= mn and x <= mx:
     1119                    j += 1
     1120            t = TimeSeries(j)
     1121            j = 0
     1122            for i from 0 <= i < self._length:
     1123                x = self._values[i]
     1124                if x >= mn and x <= mx:
     1125                    t._values[j] = x
     1126                    j += 1
     1127        return t
     1128
    10541129    def histogram(self, Py_ssize_t bins=50, bint normalize=False):
    10551130        """
    10561131        Return the frequency histogram of the values in
    cdef class TimeSeries: 
    11181193               
    11191194        INPUT:
    11201195            bins -- positive integer (default: 50)
    1121             normalize -- whether to normalize so the total area in the
    1122                          bars of the histogram is 1.
     1196            normalize -- (default: True) whether to normalize so the total
     1197                         area in the bars of the histogram is 1.
    11231198            **kwds -- passed to the bar_chart function
    11241199        OUTPUT:
    11251200            a histogram plot
    cdef class TimeSeries: 
    11531228
    11541229    def randomize(self, distribution='uniform', loc=0, scale=1, **kwds):
    11551230        """
     1231        Randomize the entries in this time series.
     1232       
    11561233        INPUT:
    1157             distribution -- 'uniform'
    1158                             'normal'
    1159                             'semicircle'
     1234            distribution -- 'uniform':    from loc to loc + scale
     1235                            'normal':     mean loc and standard deviation scale
     1236                            'semicircle': with center at loc (scale is ignored)
     1237            loc   -- float (default: 0)
     1238            scale -- float (default: 1)
     1239
     1240        NOTE: All random numbers are generated using the high quality
     1241        GMP random number function gmp_urandomb_ui.  This respects the
     1242        Sage set_random_state command.  It's not quite as fast at the
     1243        C library random number generator, but is better quality, and
     1244        is platform independent.
    11601245        """
    11611246        if distribution == 'uniform':
    11621247            self._randomize_uniform(loc, scale)
    cdef class TimeSeries: 
    11701255    def _randomize_uniform(self, double left, double right):
    11711256        if left >= right:
    11721257            raise ValueError, "left must be less than right"
     1258
     1259        cdef randstate rstate = current_randstate()
    11731260        cdef Py_ssize_t k
    11741261        cdef double d = right - left
    11751262        for k from 0 <= k < self._length:
    1176             self._values[k] = ranf() * d - left
     1263            self._values[k] = rstate.c_rand_double() * d + left
    11771264
    11781265    def _randomize_normal(self, double m, double s):
    11791266        """
    cdef class TimeSeries: 
    11851272        """
    11861273        # Ported from http://users.tkk.fi/~nbeijar/soft/terrain/source_o2/boxmuller.c
    11871274        # This the box muller algorithm.
     1275        cdef randstate rstate = current_randstate()
     1276
    11881277        cdef double x1, x2, w, y1, y2
    11891278        cdef Py_ssize_t k
    11901279        for k from 0 <= k < self._length:
    11911280            while 1:
    1192                 x1 = 2.0 * ranf() - 1.0
    1193                 x2 = 2.0 * ranf() - 1.0
     1281                x1 = 2.0 * rstate.c_rand_double() - 1.0
     1282                x2 = 2.0 * rstate.c_rand_double() - 1.0
    11941283                w = x1 * x1 + x2 * x2
    11951284                if w < 1.0: break
    11961285            w = sqrt( (-2.0 * log( w ) ) / w )
    cdef class TimeSeries: 
    12041293    def _randomize_semicircle(self, double center):
    12051294        cdef Py_ssize_t k
    12061295        cdef double x, y, s, d = 2, left = center - 1, z
     1296        cdef randstate rstate = current_randstate()
    12071297        z = d*d
    12081298        s = 1.5707963267948966192  # pi/2
    12091299        for k from 0 <= k < self._length:
    12101300            while 1:
    1211                 x = ranf() * d - 1
    1212                 y = ranf() * s
     1301                x = rstate.c_rand_double() * d - 1
     1302                y = rstate.c_rand_double() * s
    12131303                if y*y + x*x < 1:
    12141304                    break
    12151305            self._values[k] = x + center