Ticket #3356: sage-3356-part5.patch

File sage-3356-part5.patch, 44.4 KB (added by was, 14 years ago)
  • sage/finance/all.py

    # HG changeset patch
    # User William Stein <wstein@gmail.com>
    # Date 1213145628 25200
    # Node ID 3b10b5f967cbdcdca5d508fab509203a422edd69
    # Parent  3a4898af9bfb4785202cd250fa13f317d9d65054
    trac #3556 -- Finance  -- tons of additional polish, code, bug fixing, etc.
    
    diff -r 3a4898af9bfb -r 3b10b5f967cb sage/finance/all.py
    a b from stock import Stock 
    22
    33from markov_multifractal import MarkovSwitchingMultifractal
    44
    5 from time_series import TimeSeries
     5from time_series import TimeSeries, linear_filter
     6
     7from fractal import *
  • new file sage/finance/fractal.pyx

    diff -r 3a4898af9bfb -r 3b10b5f967cb sage/finance/fractal.pyx
    - +  
     1"""
     2Multifractal Random Walk
     3
     4This module implements Mandelbrot's the multifractal random walk model
     5of asset returns as developed by Bacry, Kozhemyak, and Muzy.
     6See
     7
     8"""
     9
     10from sage.rings.all import RDF, CDF, Integer
     11from sage.modules.all import vector
     12I = CDF.gen()
     13
     14from time_series cimport TimeSeries
     15
     16cdef extern from "math.h":
     17    double exp(double)
     18    double log(double)
     19    double pow(double, double)
     20    double sqrt(double)
     21
     22##################################################################
     23# Simultation
     24##################################################################
     25
     26def stationary_gaussian_simulation(s, N, n=1):
     27    """
     28    Implementation of the Davies-Harte algorithm which given an
     29    autocovariance sequence (ACVS) s and an integer N, simulates N
     30    steps of the corresponding stationary Gaussian process with mean
     31    0. We assume that a certain Fourier transform associated to s is
     32    nonnegative; if it isn't, this algorithm fails with a
     33    NotImplementedError.
     34
     35    INPUT:
     36        s -- a list of real numbers that defines the ACVS.
     37             Optimally s should have length N+1; if not
     38             we pad it with extra 0's until it has length N+1.
     39        N -- a positive integer
     40       
     41    OUTPUT:
     42        a list of n time series
     43       
     44    EXAMPLES:
     45    We define an autocovariance sequence:
     46        sage: N = 2^15
     47        sage: s = [1/math.sqrt(k+1) for k in [0..N]]
     48        sage: s[:5]
     49        [1.0, 0.70710678118654746, 0.57735026918962584, 0.5, 0.44721359549995793]
     50
     51    We run the simulation:
     52        sage: set_random_seed(0)
     53        sage: sim = finance.stationary_gaussian_simulation(s, N)[0]
     54
     55    Note that indeed the autocovariance sequence approximates s well:
     56        sage: [sim.autocovariance(i) for i in [0..4]]
     57        [0.98665816086255..., 0.69201577095377..., 0.56234006792017..., 0.48647965409871..., 0.43667043322102...]
     58
     59    WARNING: If you were to do the above computation with a small
     60    value of N, then the autocovariance sequence would not approximate
     61    s very well.
     62
     63    REFERENCES:
     64    This is a standard algorithm that is described in several papers.
     65    It is summarized nicely with many applications at the beginning of
     66    'SIMULATING A CLASS OF STATIONARY GAUSSIAN PROCESSES USING THE
     67    DAVIES-HARTE ALGORITHM, WITH APPLICATION TO LONG MEMORY
     68    PROCESSES', 2000, Peter F. Craigmile, which is easily found as a
     69    free pdf via a Google search.  This paper also generalizes the
     70    algorithm to the case when all elements of s are nonpositive.
     71
     72    The book 'Wavelet Methods for Time Series Analysis' by Percival
     73    and Walden also describes this algorithm, but has a typo in that
     74    they put a 2*pi instead of pi a certain sum.  That book describes
     75    exactly how to use Fourier transform.  The description is in
     76    Section 7.8.  Note that these pages are missing from the Google
     77    Books version of the book, but are in the Amazon.com preview of
     78    the book.
     79    """
     80    N = Integer(N)
     81    if N < 1:
     82        raise ValueError, "N must be positive"
     83
     84    if not isinstance(s, TimeSeries):
     85        s = TimeSeries(s)
     86
     87    # Make sure s has length N+1.
     88    if len(s) > N + 1:
     89        s = s[:N+1]
     90    elif len(s) < N + 1:
     91        s += TimeSeries(N+1-len(s))
     92
     93    # Make symmetrized vector.
     94    v = s + s[1:-1].reversed()
     95
     96    # Compute its fast fourier transform.
     97    cdef TimeSeries a
     98    a = v.fft()
     99
     100    # Take the real entries in the result.
     101    a = TimeSeries([a[0]]) + a[1:].scale_time(2)
     102   
     103    # Verify the nonnegativity condition.
     104    if a.min() < 0:
     105        raise NotImplementedError, "Stationary Gaussian simulation only implemented when Fourier transform is nonnegative"
     106   
     107    sims = []
     108    cdef Py_ssize_t i, k, iN = N
     109    cdef TimeSeries y, Z
     110    cdef double temp, nn = N
     111    for i from 0 <= i < n:
     112        Z = TimeSeries(2*iN).randomize('normal')
     113        y = TimeSeries(2*iN)
     114        y._values[0] = sqrt(2*nn*a._values[0]) * Z._values[0]
     115        y._values[2*iN-1] = sqrt(2*nn*a._values[iN]) * Z._values[2*iN-1]
     116        for k from 1 <= k < iN:
     117            temp = sqrt(nn*a._values[k])
     118            y._values[2*k-1] = temp*Z._values[2*k-1]
     119            y._values[2*k] = temp*Z._values[2*k]
     120        y.ifft(overwrite=True)
     121        sims.append(y[:iN])
     122
     123    return sims
     124
     125def fractional_gaussian_noise_simulation(double H, double sigma2, N, n=1):
     126    """
     127    Return $n$ simulations of with N steps each of fractional Gaussian
     128    noise with Hurst parameter $H$ and innovations variance sigma2.
     129   
     130    INPUT:
     131        H -- float; 0 < H < 1; the Hurst parameter
     132        sigma2 - float; innovation variance
     133        N -- positive integer
     134        n -- positive integer (default: 1)
     135
     136    OUTPUT:
     137        list of n time series.
     138
     139    EXAMPLES:
     140    We simulate a fractional Gaussian noise.
     141        sage: set_random_seed(0)
     142        sage: finance.fractional_gaussian_noise_simulation(0.8,1,10,2)
     143        [[-0.1157, 0.7025, 0.4949, 0.3324, 0.7110, 0.7248, -0.4048, 0.3103, -0.3465, 0.2964],
     144         [-0.5981, -0.6932, 0.5947, -0.9995, -0.7726, -0.9070, -1.3538, -1.2221, -0.0290, 1.0077]]
     145
     146    The sums define a fractional Brownian motion process:
     147        sage: set_random_seed(0)
     148        sage: finance.fractional_gaussian_noise_simulation(0.8,1,10,1)[0].sums()
     149        [-0.1157, 0.5868, 1.0818, 1.4142, 2.1252, 2.8500, 2.4452, 2.7555, 2.4090, 2.7054]
     150
     151    ALGORITHM: See 'SIMULATING A CLASS OF STATIONARY GAUSSIAN
     152    PROCESSES USING THE DAVIES-HARTE ALGORITHM, WITH APPLICATION TO
     153    LONG MEMORY PROCESSES', 2000, Peter F. Craigmile for a discussion
     154    and references for why the algorithm we give -- which uses
     155    the stationary_gaussian_simulation
     156    """
     157    if H <= 0 or H >= 1:
     158        raise ValueError, "H must satisfy 0 < H < 1"
     159    if sigma2 <= 0:
     160        raise ValueError, "sigma2 must be positive"
     161    N = Integer(N)
     162    if N < 1:
     163        raise ValueError, "N must be positive"
     164    cdef TimeSeries s = TimeSeries(N+1)
     165    s._values[0] = sigma2
     166    cdef Py_ssize_t k
     167    cdef double H2 = 2*H
     168    for k from 1 <= k <= N:
     169        s._values[k] = sigma2/2 * (pow(k+1,H2) - 2*pow(k,H2) + pow(k-1,H2))
     170    return stationary_gaussian_simulation(s, N, n)
     171
     172
     173def fractional_brownian_motion_simulation(double H, double sigma2, N, n=1):
     174    """
     175    Returns the partial sums of a fractional Gaussian noise simulation
     176    with the same input parameters.
     177
     178    INPUT:
     179        H -- float; 0 < H < 1; the Hurst parameter
     180        sigma2 - float; innovation variance
     181        N -- positive integer
     182        n -- positive integer (default: 1)
     183    OUTPUT:
     184        list of n time series.
     185
     186    EXAMPLES:
     187        sage: set_random_seed(0)
     188        sage: finance.fractional_brownian_motion_simulation(0.8,0.1,8,1)
     189        [[-0.0754, 0.2629, 0.0861, 0.2324, 0.1765, -0.0557, 0.0198, -0.0177]]
     190        sage: set_random_seed(0)
     191        sage: finance.fractional_brownian_motion_simulation(0.8,0.01,8,1)
     192        [[-0.0239, 0.0831, 0.0272, 0.0735, 0.0558, -0.0176, 0.0063, -0.0056]]
     193        sage: finance.fractional_brownian_motion_simulation(0.8,0.01,8,2)
     194        [[-0.0167, 0.0510, -0.0081, 0.0595, 0.0878, 0.0806, -0.1131, 0.0283],
     195         [0.0244, -0.0397, 0.0278, -0.0489, 0.1127, 0.0245, 0.0589, 0.0535]]
     196    """
     197    return fractional_gaussian_noise_simulation(H,sigma2,N,n)
     198   
     199def multifractal_cascade_random_walk_simulation(double T, double lambda2,
     200                                                   double ell, N, n=1):
     201    """
     202    Return a simulation of a multifractal random walk using the
     203    log-normal cascade model of Bacry-Kozhemyak-Muzy 2008.  This
     204    walk can be interpreted as the sequence of logarithms of a
     205    price series.
     206
     207    INPUT:
     208        T       -- positive real; the integral scale
     209        lambda2 -- positive real; the intermittency coefficient
     210        ell     -- a small number -- time step size.
     211        N       -- number of steps in simulation
     212        n       -- the number of simulations
     213       
     214    OUTPUT:
     215        list of time series
     216
     217    EXAMPLES:
     218        sage: set_random_seed(0)
     219        sage: a = finance.multifractal_cascade_random_walk_simulation(3770,0.02,0.01,10,3)
     220        sage: a
     221        [[-0.0096, 0.0025, 0.0066, 0.0016, 0.0078, 0.0051, 0.0047, -0.0013, 0.0003, -0.0043],
     222         [0.0003, 0.0035, 0.0257, 0.0358, 0.0377, 0.0563, 0.0661, 0.0746, 0.0749, 0.0689],
     223         [-0.0120, -0.0116, 0.0043, 0.0078, 0.0115, 0.0018, 0.0085, 0.0005, 0.0012, 0.0060]]
     224
     225    The corresponding price series:
     226        sage: a[0].exp()
     227        [0.9905, 1.0025, 1.0067, 1.0016, 1.0078, 1.0051, 1.0047, 0.9987, 1.0003, 0.9957]
     228    """
     229    if ell <= 0:
     230        raise ValueError, "ell must be positive"
     231    if T <= ell:
     232        raise ValueError, "T must be > ell"
     233    if lambda2 <= 0:
     234        raise ValueError, "lambda2 must be positive"
     235    N = Integer(N)
     236    if N < 1:
     237        raise ValueError, "N must be positive"
     238
     239    # Compute the mean of the Gaussian stationary process omega.
     240    # See page 3 of Bacry, Kozhemyak, Muzy, 2008 -- "Log-Normal
     241    # Continuous Cascades..."
     242    cdef double mean = -lambda2 * (log(T/ell) + 1)
     243
     244    # Compute the autocovariance sequence of the Gaussian
     245    # stationary process omega.  [Loc. cit.]  We have
     246    # tau = ell*i in that formula (6).
     247    cdef TimeSeries s = TimeSeries(N+1)
     248    s._values[0] = lambda2*(log(T/ell) + 1)
     249    cdef Py_ssize_t i
     250    for i from 1 <= i < N+1:
     251        if ell*i >= T:
     252            s._values[i] = lambda2 * log(T/(ell*i))
     253        else:
     254            break  # all covariance numbers after this point are 0
     255
     256    # Compute n simultations of omega, but with mean 0.
     257    omega = stationary_gaussian_simulation(s, N+1, n)
     258
     259    # Increase each by the given mean.
     260    omega = [t.add_scalar(mean) for t in omega]
     261
     262    # For each simultation, create corresponding multifractal
     263    # random walk, as explained on page 6 of [loc. cit.]
     264    sims = []
     265    cdef Py_ssize_t k
     266    cdef TimeSeries eps, steps, om
     267    for om in omega:
     268        # First compute N Gaussian white noise steps
     269        eps = TimeSeries(N).randomize('normal', 0, ell)
     270
     271        # Compute the steps of the multifractal random walk.
     272        steps = TimeSeries(N)
     273        for k from 0 <= k < N:
     274            steps._values[k] = eps._values[k] * exp(om._values[k])
     275
     276        sims.append(steps.sums())
     277       
     278    return sims
     279       
     280   
     281
     282   
     283   
     284
     285
     286## ##################################################################
     287## # Forecasting
     288## ##################################################################
     289## def multifractal_cascade_linear_filter(double T,  double lambda2,
     290##                    double sigma, double ell, Py_ssize_t M):
     291##     """
     292##     NOTE: I CANT GET ANALYTIC PARAMETERS RIGHT -- SWITCH TO MONTE CARLO...
     293   
     294##     Create the linear filter for the multifractal cascade with M steps
     295##     with given multifractal parameters.  Use this to predict the
     296##     squares of log price returns, i.e., the volatility.
     297   
     298##     INPUT:
     299##         T       -- positive real; the integral scale
     300##         lambda2 -- positive real; the intermittency coefficient
     301##         sigma   -- scaling parameter
     302##         ell     -- a small number -- time step size.
     303##         M       -- number of steps in linear filter
     304
     305##     OUTPUT:
     306##         a time series, which should be given as the input
     307##         to the linear_forecast function.
     308
     309##     REFERENCES:  See Section 13.6 "Linear Prediction and Linear Predictive Coding"
     310##     from Numerical Recipes in C for computing linear prediction
     311##     filters.  See the top of page 3 -- equation 39 -- of Bacry,
     312##     Kozhemyak, Muzy, 2006, Continuous cascade models for asset returns
     313##     for the specific analytic formulas we use of the theoretical
     314##     autocovariance.
     315
     316##     ALGORITHM: The runtime is dominated by solving a numerical linear
     317##     equation involving an M x M matrix. Thus the filter for M up to
     318##     5000 can reasonably be done in less than a minute.
     319
     320##     EXAMPLES:
     321##     We compute a linear filter with just a few terms:
     322##         sage: L = finance.multifractal_cascade_linear_filter(1000, 0.02, 2, 0.01, 5); L
     323##         [0.6453, 0.0628, 0.0822, 0.0570, 0.0959]
     324
     325##     Note that the sum as we increase the length of the filter goes to 1:
     326##         sage: sum(L)
     327##         0.94318436639725745
     328##         sage: L = finance.multifractal_cascade_linear_filter(1000, 0.02, 2, 0.01, 1000); L
     329##         [0.6014, 0.0408, 0.0576, 0.0323, 0.0241 ... 0.0001, 0.0001, 0.0002, 0.0002, 0.0005]
     330##         sage: sum(L)
     331##         0.99500051396305267
     332
     333##     We create a random multifractal walk:
     334##         sage: set_random_seed(0)
     335##         sage: y = finance.multifractal_cascade_random_walk_simulation(3770,0.02,0.01,2000,1)[0]; y
     336##         [0.0002, -0.0035, -0.0107, 0.0234, 0.0181 ... 0.2526, 0.2427, 0.2508, 0.2530, 0.2530]
     337
     338##     We then square each term in the series.
     339##         sage: y2 = y.pow(2)
     340##         sage: y2
     341##         [0.0000, 0.0000, 0.0001, 0.0005, 0.0003 ... 0.0638, 0.0589, 0.0629, 0.0640, 0.0640]
     342##         sage: y2[:-1].linear_forecast(L)
     343##         0.064481359234932076
     344##         sage: y2[:-2].linear_forecast(L)
     345##         0.063950875470110996
     346   
     347##     """
     348##     cdef TimeSeries c
     349##     cdef Py_ssize_t i
     350
     351##     if M <= 0:
     352##         raise ValueError, "M must be positive"
     353
     354##     # 1. Compute the autocovariance sequence
     355##     c = TimeSeries(M+1)
     356
     357##     for i from 1 <= i <= M:
     358##         c._values[i] = lambda2 * log(T/i)
     359##     print c
     360
     361##     # cdef double s4 = pow(sigma,4)
     362##     # cdef double ell2 = ell*ell
     363##     # cdef e = -4*lambda2
     364##     #for i from 1 <= i <= M:
     365##     #   c(i) = sigma^4 * ell^2 * (i/T)^(-4*lambda^2)
     366##     #    c._values[i] = s4 * ell2 * pow(i/T, e)
     367
     368##     #cdef double K, tau
     369##     #for i from 0 <= i <= M:
     370##     #    tau = ell*i
     371##     #    K = s4 * pow(T,4*lambda2) / ((1 + e) * (2 + e))
     372##     #    c._values[i] = K * (pow(abs(ell + tau), 2+e) + pow(abs(ell-tau), 2+e) - 2*pow(abs(tau),2+e))
     373
     374##     # Also create the numpy vector v with entries c(1), ..., c(M).
     375##     v = c[1:].numpy()
     376   
     377##     # 2. Create the autocovariances numpy 2d array A whose i,j entry
     378##     # is c(|i-j|).
     379##     import numpy
     380##     A = numpy.array([[c[abs(j-k)] for j in range(M)] for k in range(M)])
     381   
     382##     # 3. Solve the equation A * x = v
     383##     x = numpy.linalg.solve(A, v)
     384
     385##     # 4. The entries of x give the linear filter coefficients.
     386##     return TimeSeries(x)
     387
     388
     389
     390##################################################################
     391# Parameter Estimation
     392##################################################################
  • sage/finance/markov_multifractal.py

    diff -r 3a4898af9bfb -r 3b10b5f967cb sage/finance/markov_multifractal.py
    a b class MarkovSwitchingMultifractal: 
    3434            Markov switching multifractal model with m0 = 1.4, sigma = 0.5, b = 3.0, and gamma_8 = 0.95
    3535            sage: yen_usd = finance.MarkovSwitchingMultifractal(10,1.448,0.461,0.998,3.76)
    3636            sage: cad_usd = finance.MarkovSwitchingMultifractal(10,1.278,0.262,0.644,2.11)
     37            sage: dm = finance.MarkovSwitchingMultifractal(10,1.326,0.643,0.959,2.7)
    3738        """
    3839        self.__m0 = float(m0)
    3940        assert self.__m0 >= 0 and self.__m0 <= 2, "m0 must be between 0 and 2"
  • sage/finance/time_series.pyx

    diff -r 3a4898af9bfb -r 3b10b5f967cb sage/finance/time_series.pyx
    a b cdef extern from "string.h": 
    4747cdef extern from "string.h":
    4848    void* memcpy(void* dst, void* src, size_t len)
    4949
     50cdef extern from "arrayobject.h":
     51    cdef enum:
     52        NPY_OWNDATA = 0x0004 #bit mask so numpy does not free array contents when its destroyed
     53    object PyArray_FromDimsAndData(int,int*,int,double *)
     54    ctypedef int intp
     55    void import_array()
     56    ctypedef extern class numpy.ndarray [object PyArrayObject]:
     57        cdef int nd
     58        cdef int flags
     59        cdef intp *dimensions
     60        cdef char* data
     61
    5062from sage.rings.integer import Integer
    5163from sage.rings.real_double import RDF
    5264from sage.modules.real_double_vector cimport RealDoubleVectorSpaceElement
    53 
    5465
    5566max_print = 10
    5667digits = 4
    cdef class TimeSeries: 
    7889        This implicity calls init.
    7990            sage: finance.TimeSeries([pi, 3, 18.2])
    8091            [3.1416, 3.0000, 18.2000]
     92
     93        Conversion from a numpy 1d array, which is very fast.
     94            sage: v = finance.TimeSeries([1..5])
     95            sage: w = v.numpy()
     96            sage: finance.TimeSeries(w)
     97            [1.0000, 2.0000, 3.0000, 4.0000, 5.0000]
     98
     99        Conversion from an n-dimensional numpy array also works:
     100            sage: import numpy
     101            sage: v = numpy.array([[1,2], [3,4]], dtype=float); v
     102            array([[ 1.,  2.],
     103                   [ 3.,  4.]])
     104            sage: finance.TimeSeries(v)
     105            [1.0000, 2.0000, 3.0000, 4.0000]
    81106        """
    82107        cdef RealDoubleVectorSpaceElement z
     108        cdef ndarray np
    83109        if isinstance(values, (int, long, Integer)):
    84110            values = [0.0]*values
     111        elif PY_TYPE_CHECK(values, ndarray):
     112            np = values
     113            if np.nd != 1:
     114                np = values.reshape([values.size])
     115            self._length = np.dimensions[0]
     116            self._values = <double*> sage_malloc(sizeof(double) * self._length)
     117            if self._values == NULL:
     118                raise MemoryError
     119            memcpy(self._values, np.data, sizeof(double)*self._length)
     120            return
    85121        elif PY_TYPE_CHECK(values, RealDoubleVectorSpaceElement):
    86122            # Fast constructor from real double vector.
    87123            z = values
    cdef class TimeSeries: 
    465501            memcpy(v._values + i*T._length, T._values, sizeof(double)*T._length)
    466502        return v
    467503
     504    def linear_filter(self, M):
     505        acvs = [self.autocovariance(i) for i in range(M+1)]
     506        return linear_filter(acvs)
     507
     508    def linear_forecast(self, filter):
     509        cdef TimeSeries filt
     510        if isinstance(filter, TimeSeries):
     511            filt = filter
     512        else:
     513            filt = TimeSeries(filter)
     514       
     515        cdef double f = 0
     516        cdef Py_ssize_t i
     517        for i from 0 <= i < min(self._length, filt._length):
     518            f += self._values[self._length - i - 1] * filt._values[i]
     519        return f
     520
     521    def reversed(self):
     522        """
     523        Return new time series obtain from this time series by
     524        reversing the order of the entries in this time series.
     525
     526        OUTPUT:
     527            time series
     528
     529        EXAMPLES:
     530            sage: v = finance.TimeSeries([1..5])
     531            sage: v.reversed()
     532            [5.0000, 4.0000, 3.0000, 2.0000, 1.0000]
     533        """
     534        cdef Py_ssize_t i, n = self._length-1
     535        cdef TimeSeries t = new_time_series(self._length)
     536
     537        for i from 0 <= i < self._length:
     538            t._values[i] = self._values[n - i]
     539        return t       
     540
    468541    def extend(self, right):
    469542        """
    470543        Extend this time series by appending elements from the iterable right.
    cdef class TimeSeries: 
    579652            t._values[i] = self._values[i] if self._values[i] >= 0 else -self._values[i]
    580653        return t
    581654
    582     def diffs(self):
     655    def diffs(self, Py_ssize_t k = 1):
    583656        """
    584657        Return the new time series got by taking the differences of
    585658        successive terms in the time series.  So if self is the time
    586659        series $X_0, X_1, X_2, ...$, then this function outputs the
    587660        series $X_1 - X_0, X_2 - X_1, ...$.  The output series has one
    588         less term than the input series.
     661        less term than the input series.  If the optional parameter
     662        $k$ is given, return $X_k - X_0, X_{2k} - X_k, ...$.
     663
     664        INPUT:
     665            k -- positive integer (default: 1)
    589666
    590667        OUTPUT:
    591668            a new time series.
    cdef class TimeSeries: 
    596673            sage: v.diffs()
    597674            [-1.0000, -2.7000, 0.7000, 6.0000]
    598675        """
     676        if k != 1:
     677            return self.scale_time(k).diffs()
    599678        cdef Py_ssize_t i
    600679        cdef TimeSeries t = new_time_series(self._length - 1)
    601680        for i from 1 <= i < self._length:
    cdef class TimeSeries: 
    626705            sage: v.scale_time(10)
    627706            []
    628707
     708        A series of odd length:
     709            sage: v = finance.TimeSeries([1..5]); v
     710            [1.0000, 2.0000, 3.0000, 4.0000, 5.0000]
     711            sage: v.scale_time(2)
     712            [1.0000, 3.0000, 5.0000]
     713
    629714        TESTS:
    630715            sage: v.scale_time(0)
    631716            Traceback (most recent call last):
    cdef class TimeSeries: 
    639724        if k <= 0:
    640725            raise ValueError, "k must be positive"
    641726
    642         cdef Py_ssize_t i
    643         cdef TimeSeries t = new_time_series(self._length / k)  # in C / is floor division.
    644         for i from 0 <= i < self._length/k:
     727        cdef Py_ssize_t i, n
     728        n = self._length/k
     729        if self._length % 2:
     730            n += 1
     731        cdef TimeSeries t = new_time_series(n)  # in C / is floor division.
     732        for i from 0 <= i < n:
    645733            t._values[i] = self._values[i*k]
    646734        return t
    647735
    cdef class TimeSeries: 
    739827
    740828        return v
    741829
     830    def show(self, *args, **kwds):
     831        return self.plot(*args, **kwds)
    742832
    743833    def plot(self, Py_ssize_t plot_points=1000, points=False, **kwds):
    744834        """
    cdef class TimeSeries: 
    772862            s = self._length/plot_points
    773863            if plot_points <= 0:
    774864                raise ValueError, "plot_points must be a positive integer"
    775             v = self.scale_time(s).list()
     865            v = self.scale_time(s).list()[:plot_points]
    776866        else:
    777867            s = 1
    778868            v = self.list()
    cdef class TimeSeries: 
    9511041        """
    9521042        return self.sum() / self._length
    9531043
    954     def power(self, double k):
     1044    def pow(self, double k):
    9551045        """
    9561046        Return new time series with every elements of self raised to the
    9571047        kth power.
     1048
     1049        INPUT:
     1050            k -- float
     1051        OUTPUT:
     1052            time series
     1053
     1054        EXAMPLES:
     1055            sage: v = finance.TimeSeries([1,1,1,2,3]); v
     1056            [1.0000, 1.0000, 1.0000, 2.0000, 3.0000]
     1057            sage: v.pow(2)
     1058            [1.0000, 1.0000, 1.0000, 4.0000, 9.0000]       
    9581059        """
    9591060        cdef Py_ssize_t i
    9601061        cdef TimeSeries t = new_time_series(self._length)
    cdef class TimeSeries: 
    9791080            sage: v.moment(1)
    9801081            1.6000000000000001
    9811082            sage: v.moment(2)
    982             ?
     1083            3.2000000000000002
    9831084        """
    9841085        if k <= 0:
    9851086            raise ValueError, "k must be positive"
    cdef class TimeSeries: 
    10111112        # We could make this slightly faster by doing the add scalar
    10121113        # and moment calculation together.  But that would be nasty.
    10131114        return self.add_scalar(-mu).moment(k)
     1115
     1116    def covariance(self, TimeSeries other):
     1117        r"""
     1118        Return the covariance of the time series self and other.
     1119
     1120        INPUT:
     1121            self, other -- time series
     1122
     1123        Whichever time series has more terms is truncated.
     1124       
     1125        EXAMPLES:
     1126            sage: v = finance.TimeSeries([1,-2,3]); w = finance.TimeSeries([4,5,-10])
     1127            sage: v.covariance(w)
     1128            -11.777777777777779
     1129        """
     1130        cdef double mu = self.mean(), mu2 = other.mean()
     1131        cdef double s = 0
     1132        cdef Py_ssize_t i
     1133        cdef Py_ssize_t n = min(self._length, other._length)
     1134        for i from 0 <= i < n:
     1135            s += (self._values[i] - mu)*(other._values[i] - mu2)
     1136        return s / n
     1137        # NOTE: There is also a formula
     1138        #       self._values[i]*other._values[i] - mu*mu2
     1139        # but that seems less numerically stable (?), and when tested
     1140        # was not noticeably faster.
     1141
     1142    def autocovariance(self, Py_ssize_t k=0):
     1143        r"""
     1144        Return the k-th autocovariance function $\gamma(k)$ of self.
     1145        This is the covariance of self with self shifted by $k$, i.e.,
     1146        $$
     1147        ( \sum_{t=0}^{n-k-1} (x_t - \mu)(x_{t + k} - \mu) ) / n,
     1148        $$
     1149        where $n$ is the length of self.
     1150
     1151        Note the denominator of $n$, which gives a "better" sample
     1152        estimatator.
     1153
     1154        INPUT:
     1155            k -- a nonnegative integer (default: 0)
     1156
     1157        OUTPUT:
     1158            float
     1159
     1160        EXAMPLES:
     1161            sage: v = finance.TimeSeries([13,8,15,4,4,12,11,7,14,12])
     1162            sage: v.autocovariance(0)
     1163            14.4
     1164            sage: mu = v.mean(); sum([(a-mu)^2 for a in v])/len(v)
     1165            14.4
     1166            sage: v.autocovariance(1)
     1167            -2.7000000000000002
     1168            sage: mu = v.mean(); sum([(v[i]-mu)*(v[i+1]-mu) for i in range(len(v)-1)])/len(v)
     1169            -2.7000000000000002
     1170            sage: v.autocovariance(1)
     1171            -2.7000000000000002
     1172
     1173        We illustrate with a random sample that an independently and
     1174        identically distributed distribution with zero mean and
     1175        variance $\sigma^2$ has autocovariance function $\gamma(h)$
     1176        with $\gamma(0) = \sigma^2$ and $\gamma(h) = 0$ for $h\neq 0$.
     1177            sage: set_random_seed(0)
     1178            sage: v = finance.TimeSeries(10^6)
     1179            sage: v.randomize('normal', 0, 5)
     1180            [3.3835, -2.0055, 1.7882, -2.9319, -4.6827 ... -5.1868, 9.2613, 0.9274, -6.2282, -8.7652]
     1181            sage: v.autocovariance(0)
     1182            24.954106897195892
     1183            sage: v.autocovariance(1)
     1184            -0.0050839047886276651
     1185            sage: v.autocovariance(2)
     1186            0.022056325329509487
     1187            sage: v.autocovariance(3)
     1188            -0.019020003743134766
     1189        """
     1190        cdef double mu = self.mean()
     1191        cdef double s = 0
     1192        cdef Py_ssize_t i
     1193        cdef Py_ssize_t n = self._length - k
     1194        for i from 0 <= i < n:
     1195            s += (self._values[i] - mu)*(self._values[i+k] - mu)
     1196        return s / self._length
     1197
     1198    def correlation(self, TimeSeries other):
     1199        """
     1200        Return the correlation of self and other, which is the
     1201        covariance of self and other divided by the product of their
     1202        standard deviation.
     1203
     1204        INPUT:
     1205            self, other -- time series
     1206           
     1207        Whichever time series has more terms is truncated.
     1208
     1209        EXAMPLES:
     1210            sage: v = finance.TimeSeries([1,-2,3]); w = finance.TimeSeries([4,5,-10])
     1211            sage: v.correlation(w)
     1212            -0.55804160922502144
     1213            sage: v.covariance(w)/(v.standard_deviation() * w.standard_deviation())
     1214            -0.55804160922502144
     1215        """
     1216        return self.covariance(other) / (self.standard_deviation() * other.standard_deviation())
     1217
     1218    def autocorrelation(self, Py_ssize_t k=1):
     1219        r"""
     1220        Return the $k$th sample autocorrelation of this time series
     1221        $x_i$.
     1222
     1223        Let $\mu$ be the sample mean.  Then the sample autocorrelation
     1224        function is
     1225           $$
     1226              \frac{\sum_{t=0}^{n-k-1} (x_t - \mu)(x_{t+k} - \mu) }
     1227                   {\sum_{t=0}^{n-1}   (x_t - \mu)^2}
     1228           $$
     1229
     1230        Note that the variance must be nonzero or you will get a
     1231        ZeroDivisionError.
     1232
     1233        INPUT:
     1234            k -- a nonnegative integer (default: 1)
     1235       
     1236        OUTPUT:
     1237            Time series
     1238
     1239        EXAMPLE:
     1240            sage: v = finance.TimeSeries([13,8,15,4,4,12,11,7,14,12])
     1241            sage: v.autocorrelation()
     1242            -0.1875
     1243            sage: v.autocorrelation(1)
     1244            -0.1875
     1245            sage: v.autocorrelation(0)
     1246            1.0
     1247            sage: v.autocorrelation(2)
     1248            -0.20138888888888887
     1249            sage: v.autocorrelation(3)
     1250            0.18055555555555555
     1251           
     1252            sage: finance.TimeSeries([1..1000]).autocorrelation()
     1253            0.997
     1254        """
     1255        return self.autocovariance(k) / self.variance(bias=True)
    10141256
    10151257    def variance(self, bias=False):
    10161258        """
    cdef class TimeSeries: 
    12231465            sage: v = finance.TimeSeries([5,4,1.3,2,8,10,3,-5]); v
    12241466            [5.0000, 4.0000, 1.3000, 2.0000, 8.0000, 10.0000, 3.0000, -5.0000]
    12251467            sage: v.histogram(3)
    1226             ([1, 5, 2],
    1227              [(-5.0, 0.00033333333333285253),
    1228               (0.00033333333333285253, 5.0006666666666657),
    1229               (5.0006666666666657, 10.000999999999998)])
     1468            ([1, 4, 3], [(-5.0, 0.0), (0.0, 5.0), (5.0, 10.0)])
    12301469        """
    12311470        if bins <= 0:
    12321471            raise ValueError, "bins must be positive"
    cdef class TimeSeries: 
    12971536            s.ymax(max(counts))
    12981537        return s
    12991538
    1300     def numpy(self):
     1539    def numpy(self, copy=True):
    13011540        """
    1302         Return numpy 1-d array corresponding to this time series.
     1541        Return version of this time series in numpy.
    13031542       
     1543        NOTE: If copy is False, return a numpy 1d array reference to
     1544        exactly the same block of memory as this time series.  This is
     1545        very very fast, and makes it easy to quickly use all
     1546        numpy/scipy functionality on self.  However, it is dangerous
     1547        because when this time series goes out of scope and is garbage
     1548        collected, the corresponding numpy reference object will point
     1549        to garbage.
     1550
     1551        INPUT:
     1552            copy -- bool (default: True)
     1553
    13041554        OUTPUT:
    13051555            a numpy 1-d array
    13061556
    13071557        EXAMPLES:
    13081558            sage: v = finance.TimeSeries([1,-3,4.5,-2])
    1309             sage: v.numpy()
     1559            sage: w = v.numpy(copy=False); w
    13101560            array([ 1. , -3. ,  4.5, -2. ])
     1561            sage: type(w)
     1562            <type 'numpy.ndarray'>
     1563            sage: w.shape
     1564            (4,)
     1565
     1566        Notice that changing w also changes v too!
     1567            sage: w[0] = 20
     1568            sage: w
     1569            array([ 20. ,  -3. ,   4.5,  -2. ])
     1570            sage: v
     1571            [20.0000, -3.0000, 4.5000, -2.0000]
     1572
     1573        If you want a separate copy do not give the copy=False option.
     1574            sage: z = v.numpy(); z
     1575            array([ 20. ,  -3. ,   4.5,  -2. ])
     1576            sage: z[0] = -10
     1577            sage: v
     1578            [20.0000, -3.0000, 4.5000, -2.0000]
    13111579        """
    1312         import numpy
    1313         #TODO: make this faster by accessing raw memory?
    1314         return numpy.array(self.list(), dtype=float)
     1580        import_array() #This must be called before using the numpy C/api or you will get segfault
     1581        cdef int dims[1]
     1582        dims[0] = self._length
     1583        cdef ndarray n = PyArray_FromDimsAndData(1, dims, 12, self._values)
     1584        if copy:
     1585            return n.copy()
     1586        else:
     1587            return n
    13151588
    13161589    def randomize(self, distribution='uniform', loc=0, scale=1, **kwds):
    13171590        """
    1318         Randomize the entries in this time series.
     1591        Randomize the entries in this time series, and return a reference
     1592        to self.  Thus this function both changes self in place, and returns
     1593        a copy of self, for convenience.
    13191594       
    13201595        INPUT:
    13211596            distribution -- 'uniform':    from loc to loc + scale
    cdef class TimeSeries: 
    13241599            loc   -- float (default: 0)
    13251600            scale -- float (default: 1)
    13261601
    1327         NOTE: All random numbers are generated using the high quality
    1328         GMP random number function gmp_urandomb_ui.  This respects the
    1329         Sage set_random_state command.  It's not quite as fast at the
    1330         C library random number generator, but is better quality, and
    1331         is platform independent.
     1602        NOTE: All random numbers are generated using algorithms that
     1603        build on the high quality GMP random number function
     1604        gmp_urandomb_ui.  Thus this function fully respects the Sage
     1605        set_random_state command.  It's not quite as fast at the C
     1606        library random number generator, but is of course much better
     1607        quality, and is platform independent.
     1608
     1609        EXAMPLES:
     1610        We generate 5 uniform random numbers in the interval [0,1]:
     1611            sage: set_random_seed(0)
     1612            sage: finance.TimeSeries(5).randomize()
     1613            [0.8685, 0.2816, 0.0229, 0.1456, 0.7314]
     1614
     1615        We generate 5 uniform random numbers from 5 to 5+2=7:
     1616            sage: set_random_seed(0)
     1617            sage: finance.TimeSeries(10).randomize('uniform', 5, 2)
     1618            [6.7369, 5.5632, 5.0459, 5.2911, 6.4628, 5.2412, 5.2010, 5.2761, 5.5813, 5.5439]
     1619
     1620        We generate 5 normal random values with mean 0 and variance 1.
     1621            sage: set_random_seed(0)
     1622            sage: finance.TimeSeries(5).randomize('normal')
     1623            [0.6767, -0.4011, 0.3576, -0.5864, -0.9365]
     1624
     1625        We generate 10 normal random values with mean 5 and variance 2.
     1626            sage: set_random_seed(0)
     1627            sage: finance.TimeSeries(10).randomize('normal', 5, 2)
     1628            [6.3534, 4.1978, 5.7153, 3.8273, 3.1269, 2.9598, 3.7484, 6.7472, 3.8986, 4.6271]
     1629
     1630        We generate 5 values using the semicircle distribution.
     1631            sage: set_random_seed(0)
     1632            sage: finance.TimeSeries(5).randomize('semicircle')
     1633            [0.7369, -0.9541, 0.4628, -0.7990, -0.4187]
     1634
     1635        We generate 1 million normal random values and create a frequency histogram.
     1636            sage: set_random_seed(0)
     1637            sage: a = finance.TimeSeries(10^6).randomize('normal')
     1638            sage: a.histogram(10)[0]
     1639            [36, 1148, 19604, 130699, 340054, 347870, 137953, 21290, 1311, 35]
     1640
     1641        We take the above values, and compute the proportion that lie within
     1642        1, 2, 3, 5, and 6 standard deviations of the mean (0):
     1643            sage: s = a.standard_deviation()
     1644            sage: len(a.clip(-s,s))/float(len(a))
     1645            0.68309399999999998
     1646            sage: len(a.clip(-2*s,2*s))/float(len(a))
     1647            0.95455900000000005
     1648            sage: len(a.clip(-3*s,3*s))/float(len(a))
     1649            0.997228
     1650            sage: len(a.clip(-5*s,5*s))/float(len(a))
     1651            0.99999800000000005
     1652
     1653        There were no "six sigma events":
     1654            sage: len(a.clip(-6*s,6*s))/float(len(a))
     1655            1.0
    13321656        """
    13331657        if distribution == 'uniform':
    13341658            self._randomize_uniform(loc, loc + scale)
    cdef class TimeSeries: 
    13381662            self._randomize_semicircle(loc)
    13391663        else:
    13401664            raise NotImplementedError
     1665        return self
    13411666
    13421667    def _randomize_uniform(self, double left, double right):
    13431668        if left >= right:
    cdef class TimeSeries: 
    13911716                    break
    13921717            self._values[k] = x + center
    13931718
     1719
     1720    def fft(self, bint overwrite=False):
     1721        """
     1722        Return the real discrete fast fourier transform of self, as a
     1723        real time series:
     1724
     1725          [y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2))]              if n is even
     1726         
     1727          [y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2)),Im(y(n/2))]   if n is odd
     1728         
     1729        where
     1730       
     1731          y(j) = sum[k=0..n-1] x[k] * exp(-sqrt(-1)*j*k* 2*pi/n)
     1732
     1733        for j = 0..n-1.  Note that y(-j) = y(n-j).
     1734   
     1735        EXAMPLES:
     1736            sage: v = finance.TimeSeries([1..9]); v
     1737            [1.0000, 2.0000, 3.0000, 4.0000, 5.0000, 6.0000, 7.0000, 8.0000, 9.0000]
     1738            sage: w = v.fft(); w
     1739            [45.0000, -4.5000, 12.3636, -4.5000, 5.3629, -4.5000, 2.5981, -4.5000, 0.7935]
     1740
     1741        We get just the series of real parts of
     1742            sage: finance.TimeSeries([w[0]]) + w[1:].scale_time(2)
     1743            [45.0000, -4.5000, -4.5000, -4.5000, -4.5000]
     1744
     1745        An example with an even number of terms.
     1746            sage: v = finance.TimeSeries([1..10]); v
     1747            [1.0000, 2.0000, 3.0000, 4.0000, 5.0000, 6.0000, 7.0000, 8.0000, 9.0000, 10.0000]
     1748            sage: w = v.fft(); w
     1749            [55.0000, -5.0000, 15.3884, -5.0000, 6.8819, -5.0000, 3.6327, -5.0000, 1.6246, -5.0000]
     1750            sage: v.fft().ifft()
     1751            [1.0000, 2.0000, 3.0000, 4.0000, 5.0000, 6.0000, 7.0000, 8.0000, 9.0000, 10.0000]
     1752        """
     1753        import scipy.fftpack
     1754        if overwrite:
     1755            y = self
     1756        else:
     1757            y = self.__copy__()
     1758        w = y.numpy(copy=False)
     1759        scipy.fftpack.rfft(w, overwrite_x=1)
     1760        return y
     1761
     1762    def ifft(self, bint overwrite=False):
     1763        """
     1764        Return the real discrete inverse fast fourier transform of
     1765        self, which is also a real time series.
     1766
     1767        This is the inverse of fft().
     1768
     1769        The returned real array contains
     1770       
     1771                         [y(0),y(1),...,y(n-1)]
     1772                         
     1773        where for n is even
     1774       
     1775          y(j) = 1/n (sum[k=1..n/2-1] (x[2*k-1]+sqrt(-1)*x[2*k])
     1776                                       * exp(sqrt(-1)*j*k* 2*pi/n)
     1777                      + c.c. + x[0] + (-1)**(j) x[n-1])
     1778                     
     1779        and for n is odd
     1780       
     1781          y(j) = 1/n (sum[k=1..(n-1)/2] (x[2*k-1]+sqrt(-1)*x[2*k])
     1782                                       * exp(sqrt(-1)*j*k* 2*pi/n)
     1783                      + c.c. + x[0])
     1784                     
     1785        c.c. denotes complex conjugate of preceeding expression.
     1786
     1787        EXAMPLES:
     1788            sage: v = finance.TimeSeries([1..10]); v
     1789            [1.0000, 2.0000, 3.0000, 4.0000, 5.0000, 6.0000, 7.0000, 8.0000, 9.0000, 10.0000]
     1790            sage: v.ifft()
     1791            [5.1000, -5.6876, 1.4764, -1.0774, 0.4249, -0.1000, -0.2249, 0.6663, -1.2764, 1.6988]
     1792            sage: v.ifft().fft()
     1793            [1.0000, 2.0000, 3.0000, 4.0000, 5.0000, 6.0000, 7.0000, 8.0000, 9.0000, 10.0000]
     1794        """
     1795        import scipy.fftpack
     1796        if overwrite:
     1797            y = self
     1798        else:
     1799            y = self.__copy__()
     1800        w = y.numpy(copy=False)
     1801        scipy.fftpack.irfft(w, overwrite_x=1)
     1802        return y
     1803       
     1804
    13941805cdef new_time_series(Py_ssize_t length):
    13951806    """
    13961807    Return a new uninitialized time series of the given length.
    def unpickle_time_series_v1(v, Py_ssize_ 
    14381849    memcpy(t._values, PyString_AsString(v), n*sizeof(double))
    14391850    return t
    14401851
     1852
     1853
     1854
     1855def linear_filter(acvs):
     1856    """
     1857    Create a linear filter with given autococvariance sequence.
     1858   
     1859    EXAMPLES:
     1860
     1861    In this example we consider the multifractal cascade random walk
     1862    of length 1000, and use simultations to estimate the
     1863    expected first few autocovariance parameters for this model, then
     1864    use them to construct a linear filter that works vastly better
     1865    than a linear filter constructed from the same data but not using
     1866    this model. The Monte-Carlo method illustrated below should work for
     1867    predicting one "time step" into the future for any
     1868    model that can be simultated.  To predict k time steps into the
     1869    future would require using a similar technique but would require
     1870    scaling time by k.
     1871
     1872    We create 100 simultations ofa multifractal random walk.  This
     1873    models the logarithms of a stock price sequence.
     1874        sage: set_random_seed(0)
     1875        sage: y = finance.multifractal_cascade_random_walk_simulation(3700,0.02,0.01,1000,100)
     1876
     1877    For each walk below we replace the walk by the walk but where each
     1878    step size is replaced by its absolute value -- this is what we
     1879    expect to be able to predict given the model, which is only a
     1880    model for predicting volatility.  We compute the first 200
     1881    autocovariance values for every random walk:
     1882        sage: c = [[a.diffs().abs().sums().autocovariance(i) for a in y] for i in range(200)]
     1883
     1884    We make a time series out of the expected values of the
     1885    autocovariances:
     1886        sage: ac = finance.TimeSeries([finance.TimeSeries(z).mean() for z in c])
     1887        sage: ac
     1888        [3.9962, 3.9842, 3.9722, 3.9601, 3.9481 ... 1.7144, 1.7033, 1.6922, 1.6812, 1.6701]
     1889
     1890    Note: ac looks like a line -- one could best fit it to yield a lot
     1891    more approximate autocovariances.
     1892
     1893    We compute the linear filter given by the above autocovariances:
     1894        sage: F = finance.linear_filter(ac); F
     1895        [0.9982, -0.0002, -0.0002, 0.0003, 0.0001 ... 0.0002, -0.0002, -0.0000, -0.0002, -0.0014]
     1896
     1897    Note that the sum is close to 1.
     1898        sage: sum(F)
     1899        0.99593284089454...
     1900
     1901    Now we make up an 'out of sample' sequence:
     1902        sage: y2 = finance.multifractal_cascade_random_walk_simulation(3700,0.02,0.01,1000,1)[0].diffs().abs().sums()
     1903        sage: y2
     1904        [0.0013, 0.0059, 0.0066, 0.0068, 0.0184 ... 6.8004, 6.8009, 6.8063, 6.8090, 6.8339]
     1905
     1906    And we forecast the very last value using our linear filter; the forecast is close:
     1907        sage: y2[:-1].linear_forecast(F)
     1908        6.7836741372407...
     1909
     1910    In fact it is closer than we would get by forecasting using a
     1911    linear filter made from all the autocovariances of our sequence:
     1912        sage: y2[:-1].linear_forecast(y2[:-1].linear_filter(len(y2)))
     1913        6.7701687056683...
     1914
     1915    We record the last 20 forecasts, always using all correct values up to the
     1916    one we are forecasting:
     1917        sage: s1 = sum([(y2[:-i].linear_forecast(F)-y2[-i])^2 for i in range(1,20)])
     1918
     1919    We do the same, but using the autocovariances of the sample sequence:
     1920        sage: F2 = y2[:-100].linear_filter(len(F))
     1921        sage: s2 = sum([(y2[:-i].linear_forecast(F2)-y2[-i])^2 for i in range(1,20)])
     1922
     1923    Our model gives us something that is 15 percent better in this case:
     1924        sage: s2/s1
     1925        1.154646361026...
     1926
     1927    How does it compare overall?  To find out we do 100 simulations
     1928    and for each we compute the percent that our model beats naively
     1929    using the autocovariances of the sample:
     1930        sage: y_out = finance.multifractal_cascade_random_walk_simulation(3700,0.02,0.01,1000,100)
     1931        sage: s1 = []; s2 = []
     1932        sage: for v in y_out:
     1933        ...       s1.append(sum([(v[:-i].linear_forecast(F)-v[-i])^2 for i in range(1,20)]))
     1934        ...       F2 = v[:-len(F)].linear_filter(len(F))
     1935        ...       s2.append(sum([(v[:-i].linear_forecast(F2)-v[-i])^2 for i in range(1,20)]))
     1936        ...
     1937
     1938    We find that overall the model beats naive linear forecasting by 35 percent!
     1939        sage: s = finance.TimeSeries([s2[i]/s1[i] for i in range(len(s1))])
     1940        sage: s.mean()
     1941        1.354073591877...
     1942    """
     1943    cdef TimeSeries c
     1944    cdef Py_ssize_t i
     1945
     1946    M = len(acvs)-1
     1947
     1948    if M <= 0:
     1949        raise ValueError, "M must be positive"
     1950
     1951    if not isinstance(acvs, TimeSeries):
     1952        c = TimeSeries(acvs)
     1953    else:
     1954        c = acvs
     1955
     1956    # Also create the numpy vector v with entries c(1), ..., c(M).
     1957    v = c[1:].numpy()
     1958
     1959    # 2. Create the autocovariances numpy 2d array A whose i,j entry
     1960    # is c(|i-j|).
     1961    import numpy
     1962    A = numpy.array([[c[abs(j-k)] for j in range(M)] for k in range(M)])
     1963
     1964    # 3. Solve the equation A * x = v
     1965    x = numpy.linalg.solve(A, v)
     1966
     1967    # 4. The entries of x give the linear filter coefficients.
     1968    return TimeSeries(x)
  • setup.py

    diff -r 3a4898af9bfb -r 3b10b5f967cb setup.py
    a b symmetrica = Extension('sage.libs.symmet 
    475475                       include_dirs=debian_include_dirs + ['/usr/include/malloc/'],
    476476                       libraries = ["symmetrica"])
    477477
    478 time_series = Extension('sage.finance.time_series',['sage/finance/time_series.pyx'])
     478time_series = Extension('sage.finance.time_series',['sage/finance/time_series.pyx'],
     479                        include_dirs=debian_include_dirs + [SAGE_ROOT+'/local/lib/python2.5/site-packages/numpy/core/include/numpy'])
    479480
    480481markov_multifractal = Extension('sage.finance.markov_multifractal_cython',
    481482                                ['sage/finance/markov_multifractal_cython.pyx'])
    482483
     484finance_fractal = Extension('sage.finance.fractal', ['sage/finance/fractal.pyx'])
    483485
    484486#####################################################
    485487
    ext_modules = [ \ 
    596598    time_series,
    597599   
    598600    markov_multifractal,
     601    finance_fractal,
    599602
    600603    Extension('sage.media.channels',
    601604              sources = ['sage/media/channels.pyx']), \