Ticket #3356: sage-3356-part3.patch

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

    # HG changeset patch
    # User William Stein <wstein@gmail.com>
    # Date 1212541198 25200
    # Node ID a64e6d49cfac5001b3dfd7d9fb2ceab324f4db86
    # Parent  1f137599f3e3141f261529ab6930ae9880b2443d
    Fix a few tiny bugs, a memleak, and implements moments and pow method.
    
    diff -r 1f137599f3e3 -r a64e6d49cfac sage/finance/markov_multifractal.py
    a b class MarkovSwitchingMultifractal: 
    192192            sage: cad_usd = finance.MarkovSwitchingMultifractal(10,1.278,0.262,0.644,2.11); cad_usd
    193193            Markov switching multifractal model with m0 = 1.278, sigma = 0.262, b = 2.11, and gamma_10 = 0.644
    194194        """
    195         return markov_multifractal_cython.simulation(n, k,
     195        return markov_multifractal_cython.simulations(n, k,
    196196                   self.__m0, self.__sigma, self.__b,
    197197                   self.__kbar, self.gamma())
    198198       
  • sage/finance/markov_multifractal_cython.pyx

    diff -r 1f137599f3e3 -r a64e6d49cfac sage/finance/markov_multifractal_cython.pyx
    a b cdef extern from "math.h": 
    1111
    1212from time_series cimport TimeSeries
    1313
    14 def simulation(Py_ssize_t n, Py_ssize_t k,
     14def simulations(Py_ssize_t n, Py_ssize_t k,
    1515               double m0, double sigma, double b,
    1616               int kbar, gamma):
    1717    cdef double m1 = 2 - m0
    def simulation(Py_ssize_t n, Py_ssize_t  
    2929    for i from 0 <= i < k:
    3030        # Initalize the model
    3131        for j from 0 <= j < kbar:
    32             markov_state_vector._values[j] = m0 if (rstate.c_random() & 1) else m1    # n & 1 means "is odd"
     32            # n & 1 means "is odd"
     33            markov_state_vector._values[j] = m0 if (rstate.c_random() & 1) else m1   
    3334        t = TimeSeries(n)
    3435       
    3536        # Generate n normally distributed random numbers with mean 0
    def simulation(Py_ssize_t n, Py_ssize_t  
    4243            t._values[a] = sigma * eps._values[a] * sqrt(markov_state_vector.prod())
    4344
    4445            # Now update the volatility state vector
    45             j = a * kbar
    4646            for c from 0 <= c < kbar:
    4747                if rstate.c_rand_double() <= gamma_vals._values[c]:
    48                     markov_state_vector._values[k] = m0 if (rstate.c_random() & 1) else m1
     48                    markov_state_vector._values[c] = m0 if (rstate.c_random() & 1) else m1
    4949
    5050        S.append(t)
    5151       
  • sage/finance/time_series.pyx

    diff -r 1f137599f3e3 -r a64e6d49cfac sage/finance/time_series.pyx
    a b cdef extern from "math.h": 
    4141    double exp(double)
    4242    double floor(double)
    4343    double log(double)
     44    double pow(double, double)
    4445    double sqrt(double)
    4546
    4647cdef extern from "string.h":
    cdef class TimeSeries: 
    393394        the add_entries method.
    394395       
    395396        INPUT:
    396             right -- iterable that can be converted to a time series       
     397            right -- a time series
    397398        OUTPUT:
    398399            a time series
    399400
    cdef class TimeSeries: 
    401402            sage: v = finance.TimeSeries([1,2,3]); w = finance.TimeSeries([1,2])
    402403            sage: v + w
    403404            [1.0000, 2.0000, 3.0000, 1.0000, 2.0000]
    404             sage: v + xrange(4)
    405             [1.0000, 2.0000, 3.0000, 0.0000, 1.0000, 2.0000, 3.0000]
    406405            sage: v = finance.TimeSeries([1,2,-5]); v
    407406            [1.0000, 2.0000, -5.0000]
     407
     408        Note that both summands must be a time series:
     409            sage: v + xrange(4)
     410            Traceback (most recent call last):
     411            ...
     412            TypeError: right operand must be a time series
    408413            sage: [1,5] + v
    409             [1.0000, 5.0000, 1.0000, 2.0000, -5.0000]
     414            Traceback (most recent call last):
     415            ...
     416            TypeError: left operand must be a time series
    410417        """
    411418        if not isinstance(right, TimeSeries):
    412             raise TypeError, "right must be a time series"
     419            raise TypeError, "right operand must be a time series"
    413420        if not isinstance(left, TimeSeries):
    414             raise TypeError, "right must be a time series"
     421            raise TypeError, "left operand must be a time series"
    415422        cdef TimeSeries R = right
    416423        cdef TimeSeries L = left
    417424        cdef TimeSeries t = new_time_series(L._length + R._length)
    cdef class TimeSeries: 
    943950            1.6000000000000001
    944951        """
    945952        return self.sum() / self._length
     953
     954    def power(self, double k):
     955        """
     956        Return new time series with every elements of self raised to the
     957        kth power.
     958        """
     959        cdef Py_ssize_t i
     960        cdef TimeSeries t = new_time_series(self._length)
     961        for i from 0 <= i < self._length:
     962            t._values[i] = pow(self._values[i], k)
     963        return t
     964
     965    def moment(self, int k):
     966        """
     967        Return the k-th moment of self, which is just the
     968        mean of the k-th powers of the elements of self.
     969
     970        INPUT:
     971            k -- a positive integer
     972
     973        OUTPUT:
     974            double
     975           
     976        EXAMPLES:
     977            sage: v = finance.TimeSeries([1,1,1,2,3]); v
     978            [1.0000, 1.0000, 1.0000, 2.0000, 3.0000]
     979            sage: v.moment(1)
     980            1.6000000000000001
     981            sage: v.moment(2)
     982            ?
     983        """
     984        if k <= 0:
     985            raise ValueError, "k must be positive"
     986        if k == 1:
     987            return self.mean()
     988        cdef double s = 0
     989        cdef Py_ssize_t i
     990        for i from 0 <= i < self._length:
     991            s += pow(self._values[i], k)
     992        return s / self._length
     993
     994    def central_moment(self, int k):
     995        """
     996        Return the k-th central moment of self, which is just the mean
     997        of the k-th powers of the differences self[i]-mu, where mu is
     998        the mean of self.
     999
     1000        INPUT:
     1001            k -- a positive integer
     1002        OUTPUT:
     1003            double
     1004
     1005        EXAMPLES:
     1006           
     1007        """
     1008        if k == 1:
     1009            return float(0)
     1010        mu = self.mean()
     1011        # We could make this slightly faster by doing the add scalar
     1012        # and moment calculation together.  But that would be nasty.
     1013        return self.add_scalar(-mu).moment(k)
    9461014
    9471015    def variance(self, bias=False):
    9481016        """