Ticket #3356: sage-3356-part3.patch

sage/finance/markov_multifractal.py

@@ -192 +192 @@
             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
         """
-        return markov_multifractal_cython.simulation(n, k,
+        return markov_multifractal_cython.simulations(n, k,
                    self.__m0, self.__sigma, self.__b,
                    self.__kbar, self.gamma())
         

sage/finance/markov_multifractal_cython.pyx

@@ -11 +11 @@
 
 from time_series cimport TimeSeries
 
-def simulation(Py_ssize_t n, Py_ssize_t k,
+def simulations(Py_ssize_t n, Py_ssize_t k,
                double m0, double sigma, double b,
                int kbar, gamma):
     cdef double m1 = 2 - m0
@@ -29 +29 @@
     for i from 0 <= i < k:
         # Initalize the model
         for j from 0 <= j < kbar:
-            markov_state_vector._values[j] = m0 if (rstate.c_random() & 1) else m1    # n & 1 means "is odd"
+            # 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
@@ -42 +43 @@
             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 
+                    markov_state_vector._values[c] = m0 if (rstate.c_random() & 1) else m1 
 
         S.append(t)
         

sage/finance/time_series.pyx

@@ -41 +41 @@
     double exp(double)
     double floor(double)
     double log(double)
+    double pow(double, double)
     double sqrt(double)
 
 cdef extern from "string.h":
@@ -393 +394 @@
         the add_entries method.
         
         INPUT:
-            right -- iterable that can be converted to a time series        
+            right -- a time series
         OUTPUT:
             a time series
 
@@ -401 +402 @@
             sage: v = finance.TimeSeries([1,2,3]); w = finance.TimeSeries([1,2])
             sage: v + w
             [1.0000, 2.0000, 3.0000, 1.0000, 2.0000]
-            sage: v + xrange(4)
-            [1.0000, 2.0000, 3.0000, 0.0000, 1.0000, 2.0000, 3.0000]
             sage: v = finance.TimeSeries([1,2,-5]); v
             [1.0000, 2.0000, -5.0000]
+
+        Note that both summands must be a time series: 
+            sage: v + xrange(4)
+            Traceback (most recent call last):
+            ...
+            TypeError: right operand must be a time series
             sage: [1,5] + v
-            [1.0000, 5.0000, 1.0000, 2.0000, -5.0000]
+            Traceback (most recent call last):
+            ...
+            TypeError: left operand must be a time series
         """
         if not isinstance(right, TimeSeries):
-            raise TypeError, "right must be a time series"
+            raise TypeError, "right operand must be a time series"
         if not isinstance(left, TimeSeries):
-            raise TypeError, "right must be a time series"
+            raise TypeError, "left operand must be a time series"
         cdef TimeSeries R = right
         cdef TimeSeries L = left
         cdef TimeSeries t = new_time_series(L._length + R._length)
@@ -943 +950 @@
             1.6000000000000001
         """
         return self.sum() / self._length
+
+    def power(self, double k):
+        """
+        Return new time series with every elements of self raised to the
+        kth power.
+        """
+        cdef Py_ssize_t i
+        cdef TimeSeries t = new_time_series(self._length)
+        for i from 0 <= i < self._length:
+            t._values[i] = pow(self._values[i], k)
+        return t
+
+    def moment(self, int k):
+        """
+        Return the k-th moment of self, which is just the
+        mean of the k-th powers of the elements of self. 
+
+        INPUT:
+            k -- a positive integer
+
+        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.moment(1)
+            1.6000000000000001
+            sage: v.moment(2)
+            ?
+        """
+        if k <= 0:
+            raise ValueError, "k must be positive"
+        if k == 1:
+            return self.mean()
+        cdef double s = 0
+        cdef Py_ssize_t i
+        for i from 0 <= i < self._length:
+            s += pow(self._values[i], k)
+        return s / self._length
+
+    def central_moment(self, int k):
+        """
+        Return the k-th central moment of self, which is just the mean
+        of the k-th powers of the differences self[i]-mu, where mu is
+        the mean of self.
+
+        INPUT:
+            k -- a positive integer
+        OUTPUT:
+            double
+
+        EXAMPLES:
+            
+        """
+        if k == 1:
+            return float(0)
+        mu = self.mean()
+        # We could make this slightly faster by doing the add scalar
+        # and moment calculation together.  But that would be nasty.
+        return self.add_scalar(-mu).moment(k)
 
     def variance(self, bias=False):
         """