Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
""" Time Series
This is a module for working with discrete floating point time series. It is designed so that every operation is very fast, typically much faster than with other generic code, e.g., Python lists of doubles or even NumPy arrays. The semantics of time series is more similar to Python lists of doubles than Sage real double vectors or NumPy 1-D arrays. In particular, time series are not endowed with much algebraic structure and are always mutable.
.. NOTE::
NumPy arrays are faster at slicing, since slices return references, and NumPy arrays have strides. However, this speed at slicing makes NumPy slower at certain other operations.
EXAMPLES::
sage: set_random_seed(1) sage: t = finance.TimeSeries([random()-0.5 for _ in range(10)]); t [0.3294, 0.0959, -0.0706, -0.4646, 0.4311, 0.2275, -0.3840, -0.3528, -0.4119, -0.2933] sage: t.sums() [0.3294, 0.4253, 0.3547, -0.1099, 0.3212, 0.5487, 0.1647, -0.1882, -0.6001, -0.8933] sage: t.exponential_moving_average(0.7) [0.0000, 0.3294, 0.1660, 0.0003, -0.3251, 0.2042, 0.2205, -0.2027, -0.3078, -0.3807] sage: t.standard_deviation() 0.33729638212891383 sage: t.mean() -0.08933425506929439 sage: t.variance() 0.1137688493972542...
AUTHOR:
- William Stein """
#***************************************************************************** # Copyright (C) 2008 William Stein <wstein@gmail.com> # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 2 of the License, or # (at your option) any later version. # http://www.gnu.org/licenses/ #***************************************************************************** from __future__ import absolute_import
from cpython.bytes cimport PyBytes_FromStringAndSize, PyBytes_AsString from libc.math cimport exp, floor, log, pow, sqrt from libc.string cimport memcpy from cysignals.memory cimport sig_malloc, sig_free from sage.structure.richcmp cimport rich_to_bool
cimport numpy as cnumpy
from sage.misc.randstate cimport randstate, current_randstate from sage.modules.vector_real_double_dense cimport Vector_real_double_dense
cdef class TimeSeries: def __cinit__(self): """ Create new empty uninitialized time series.
EXAMPLES:
This implicitly calls new::
sage: finance.TimeSeries([1,3,-4,5]) [1.0000, 3.0000, -4.0000, 5.0000] """
def __init__(self, values, bint initialize=True): """ Initialize new time series.
INPUT:
- ``values`` -- integer (number of values) or an iterable of floats.
- ``initialize`` -- bool (default: ``True``); if ``False``, do not bother to zero out the entries of the new time series. For large series that you are going to just fill in, this can be way faster.
EXAMPLES:
This implicitly calls init::
sage: finance.TimeSeries([pi, 3, 18.2]) [3.1416, 3.0000, 18.2000]
Conversion from a NumPy 1-D array, which is very fast::
sage: v = finance.TimeSeries([1..5]) sage: w = v.numpy() sage: finance.TimeSeries(w) [1.0000, 2.0000, 3.0000, 4.0000, 5.0000]
Conversion from an n-dimensional NumPy array also works::
sage: import numpy sage: v = numpy.array([[1,2], [3,4]], dtype=float); v array([[ 1., 2.], [ 3., 4.]]) sage: finance.TimeSeries(v) [1.0000, 2.0000, 3.0000, 4.0000] sage: finance.TimeSeries(v[:,0]) [1.0000, 3.0000] sage: u = numpy.array([[1,2],[3,4]]) sage: finance.TimeSeries(u) [1.0000, 2.0000, 3.0000, 4.0000]
For speed purposes we don't initialize (so value is garbage)::
sage: t = finance.TimeSeries(10, initialize=False) """ cdef Vector_real_double_dense z cdef cnumpy.ndarray np cdef double *np_data cdef unsigned int j np = values._vector_numpy else:
# Make the array be the correct type and have a C array # for a data structure. If the array already is the # correct type and has a C array, nothing is done, so this # should be fast in the common case. raise MemoryError
else:
raise MemoryError cdef Py_ssize_t i else:
def __reduce__(self): """ Used in pickling time series.
EXAMPLES::
sage: v = finance.TimeSeries([1,-3.5]) sage: v.__reduce__() (<built-in function unpickle_time_series_v1>, (..., 2)) sage: loads(dumps(v)) == v True
Note that dumping and loading with compress ``False`` is much faster, though dumping with compress ``True`` can save a lot of space::
sage: v = finance.TimeSeries([1..10^5]) sage: loads(dumps(v, compress=False),compress=False) == v True """
def __richcmp__(TimeSeries self, other, int op): """ Compare ``self`` and ``other``. This has the same semantics as list comparison.
EXAMPLES::
sage: v = finance.TimeSeries([1,2,3]); w = finance.TimeSeries([1,2]) sage: v < w False sage: w < v True sage: v == v True sage: w == w True """ cdef Py_ssize_t i cdef double d return NotImplemented
def __dealloc__(self): """ Free up memory used by a time series.
EXAMPLES:
This tests ``__dealloc__`` implicitly::
sage: v = finance.TimeSeries([1,3,-4,5]) sage: del v """
def vector(self): """ Return real double vector whose entries are the values of this time series. This is useful since vectors have standard algebraic structure and play well with matrices.
OUTPUT:
A real double vector.
EXAMPLES::
sage: v = finance.TimeSeries([1..10]) sage: v.vector() (1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0) """ # A copy of the numpy array is made in the vector constructor
def __repr__(self): """ Return string representation of ``self``.
EXAMPLES::
sage: v = finance.TimeSeries([1,3.1908439,-4,5.93932]) sage: v.__repr__() '[1.0000, 3.1908, -4.0000, 5.9393]'
By default 4 digits after the decimal point are displayed. To change this, change ``self.finance.time_series.digits``. ::
sage: sage.finance.time_series.digits = 2 sage: v.__repr__() '[1.00, 3.19, -4.00, 5.94]' sage: v [1.00, 3.19, -4.00, 5.94] sage: sage.finance.time_series.digits = 4 sage: v [1.0000, 3.1908, -4.0000, 5.9393] """
def _repr(self, prec=None): """ Print representation of a time series.
INPUT:
- ``prec`` -- (default: ``None``) number of digits of precision or ``None``. If ``None`` use the default ``sage.finance.time_series.digits``.
OUTPUT:
A string.
EXAMPLES::
sage: v = finance.TimeSeries([1,3.1908439,-4,5.93932]) sage: v._repr() '[1.0000, 3.1908, -4.0000, 5.9393]' sage: v._repr(10) '[1.0000000000, 3.1908439000, -4.0000000000, 5.9393200000]' sage: v._repr(2) '[1.00, 3.19, -4.00, 5.94]' """ else:
def __len__(self): """ Return the number of entries in this time series.
OUTPUT:
Python integer.
EXAMPLES::
sage: v = finance.TimeSeries([1,3.1908439,-4,5.93932]) sage: v.__len__() 4 sage: len(v) 4 """
def __getitem__(self, i): """ Return i-th entry or slice of ``self``.
EXAMPLES::
sage: v = finance.TimeSeries([1,-4,3,-2.5,-4,3]) sage: v[2] 3.0 sage: v[-1] 3.0 sage: v[-10] Traceback (most recent call last): ... IndexError: TimeSeries index out of range sage: v[5] 3.0 sage: v[6] Traceback (most recent call last): ... IndexError: TimeSeries index out of range
Some slice examples::
sage: v[-3:] [-2.5000, -4.0000, 3.0000] sage: v[-3:-1] [-2.5000, -4.0000] sage: v[::2] [1.0000, 3.0000, -4.0000] sage: v[3:20] [-2.5000, -4.0000, 3.0000] sage: v[3:2] []
Make a copy::
sage: v[:] [1.0000, -4.0000, 3.0000, -2.5000, -4.0000, 3.0000]
Reverse the time series::
sage: v[::-1] [3.0000, -4.0000, -2.5000, 3.0000, -4.0000, 1.0000] """ cdef Py_ssize_t start, stop, step, j cdef TimeSeries t start = self._length - 1 else: else:
def __setitem__(self, Py_ssize_t i, double x): """ Set the i-th entry of ``self`` to ``x``.
INPUT:
- i -- a nonnegative integer.
- x -- a float.
EXAMPLES::
sage: v = finance.TimeSeries([1,3,-4,5.93932]); v [1.0000, 3.0000, -4.0000, 5.9393] sage: v[0] = -5.5; v [-5.5000, 3.0000, -4.0000, 5.9393] sage: v[-1] = 3.2; v [-5.5000, 3.0000, -4.0000, 3.2000] sage: v[10] Traceback (most recent call last): ... IndexError: TimeSeries index out of range sage: v[-5] Traceback (most recent call last): ... IndexError: TimeSeries index out of range """ raise IndexError("TimeSeries index out of range") raise IndexError("TimeSeries index out of range")
def __copy__(self): """ Return a copy of ``self``.
EXAMPLES::
sage: v = finance.TimeSeries([1,-4,3,-2.5,-4,3]) sage: v.__copy__() [1.0000, -4.0000, 3.0000, -2.5000, -4.0000, 3.0000] sage: copy(v) [1.0000, -4.0000, 3.0000, -2.5000, -4.0000, 3.0000] sage: copy(v) is v False """ cdef Py_ssize_t i
def __add__(left, right): """ Concatenate the time series ``self`` and ``right``.
.. NOTE::
To add a single number to the entries of a time series, use the ``add_scalar`` method, and to add componentwise use the ``add_entries`` method.
INPUT:
- ``right`` -- a time series.
OUTPUT:
A time series.
EXAMPLES::
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 = finance.TimeSeries([1,2,-5]); v [1.0000, 2.0000, -5.0000]
Note that both summands must be a time series::
sage: v + list(range(4)) Traceback (most recent call last): ... TypeError: right operand must be a time series sage: [1,5] + v Traceback (most recent call last): ... TypeError: left operand must be a time series """
def __mul__(left, right): """ Multiply a time series by an integer n, which (like for lists) results in the time series concatenated with itself n times.
.. NOTE::
To multiply all the entries of a time series by a single scalar, use the ``scale`` method.
INPUT:
- ``left``, ``right`` -- an integer and a time series.
OUTPUT:
A time series.
EXAMPLES::
sage: v = finance.TimeSeries([1,2,-5]); v [1.0000, 2.0000, -5.0000] sage: v*3 [1.0000, 2.0000, -5.0000, 1.0000, 2.0000, -5.0000, 1.0000, 2.0000, -5.0000] sage: 3*v [1.0000, 2.0000, -5.0000, 1.0000, 2.0000, -5.0000, 1.0000, 2.0000, -5.0000] sage: v*v # py2 Traceback (most recent call last): ... TypeError: 'sage.finance.time_series.TimeSeries' object cannot be interpreted as an index sage: v*v # py3 Traceback (most recent call last): ... TypeError: 'sage.finance.time_series.TimeSeries' object cannot be interpreted as an integer """ cdef Py_ssize_t n, i cdef TimeSeries T else: # Make n copies of T concatenated together
def autoregressive_fit(self,M): r""" This method fits the time series to an autoregressive process of order ``M``. That is, we assume the process is given by `X_t-\mu=a_1(X_{t-1}-\mu)+a_2(X_{t-1}-\mu)+\cdots+a_M(X_{t-M}-\mu)+Z_t` where `\mu` is the mean of the process and `Z_t` is noise. This method returns estimates for `a_1,\dots,a_M`.
The method works by solving the Yule-Walker equations `\Gamma a =\gamma`, where `\gamma=(\gamma(1),\dots,\gamma(M))`, `a=(a_1,\dots,a_M)` with `\gamma(i)` the autocovariance of lag `i` and `\Gamma_{ij}=\gamma(i-j)`.
.. WARNING::
The input sequence is assumed to be stationary, which means that the autocovariance `\langle y_j y_k \rangle` depends only on the difference `|j-k|`.
INPUT:
- ``M`` -- an integer.
OUTPUT:
A time series -- the coefficients of the autoregressive process.
EXAMPLES::
sage: set_random_seed(0) sage: v = finance.TimeSeries(10^4).randomize('normal').sums() sage: F = v.autoregressive_fit(100) sage: v [0.6767, 0.2756, 0.6332, 0.0469, -0.8897 ... 87.6759, 87.6825, 87.4120, 87.6639, 86.3194] sage: v.autoregressive_forecast(F) 86.0177285042... sage: F [1.0148, -0.0029, -0.0105, 0.0067, -0.0232 ... -0.0106, -0.0068, 0.0085, -0.0131, 0.0092]
sage: set_random_seed(0) sage: t=finance.TimeSeries(2000) sage: z=finance.TimeSeries(2000) sage: z.randomize('normal',1) [1.6767, 0.5989, 1.3576, 0.4136, 0.0635 ... 1.0057, -1.1467, 1.2809, 1.5705, 1.1095] sage: t[0]=1 sage: t[1]=2 sage: for i in range(2,2000): ....: t[i]=t[i-1]-0.5*t[i-2]+z[i] sage: c=t[0:-1].autoregressive_fit(2) #recovers recurrence relation sage: c #should be close to [1,-0.5] [1.0371, -0.5199] """
def autoregressive_forecast(self, filter): """ Given the autoregression coefficients as outputted by the ``autoregressive_fit`` command, compute the forecast for the next term in the series.
INPUT:
- ``filter`` -- a time series outputted by the ``autoregressive_fit`` command.
EXAMPLES::
sage: set_random_seed(0) sage: v = finance.TimeSeries(100).randomize('normal').sums() sage: F = v[:-1].autoregressive_fit(5); F [1.0019, -0.0524, -0.0643, 0.1323, -0.0539] sage: v.autoregressive_forecast(F) 11.7820298611... sage: v [0.6767, 0.2756, 0.6332, 0.0469, -0.8897 ... 9.2447, 9.6709, 10.4037, 10.4836, 12.1960] """ cdef TimeSeries filt else: filt = TimeSeries(filter)
cdef Py_ssize_t i
def reversed(self): """ Return new time series obtain from this time series by reversing the order of the entries in this time series.
OUTPUT:
A time series.
EXAMPLES::
sage: v = finance.TimeSeries([1..5]) sage: v.reversed() [5.0000, 4.0000, 3.0000, 2.0000, 1.0000] """
def extend(self, right): """ Extend this time series by appending elements from the iterable ``right``.
INPUT:
- ``right`` -- iterable that can be converted to a time series.
EXAMPLES::
sage: v = finance.TimeSeries([1,2,-5]); v [1.0000, 2.0000, -5.0000] sage: v.extend([-3.5, 2]) sage: v [1.0000, 2.0000, -5.0000, -3.5000, 2.0000] """ return raise MemoryError
def list(self): """ Return list of elements of ``self``.
EXAMPLES::
sage: v = finance.TimeSeries([1,-4,3,-2.5,-4,3]) sage: v.list() [1.0, -4.0, 3.0, -2.5, -4.0, 3.0] """ cdef Py_ssize_t i
def log(self): """ Return new time series got by taking the logarithms of all the terms in the time series.
OUTPUT:
A new time series.
EXAMPLES:
We exponentiate then log a time series and get back the original series::
sage: v = finance.TimeSeries([1,-4,3,-2.5,-4,3]); v [1.0000, -4.0000, 3.0000, -2.5000, -4.0000, 3.0000] sage: v.exp() [2.7183, 0.0183, 20.0855, 0.0821, 0.0183, 20.0855] sage: v.exp().log() [1.0000, -4.0000, 3.0000, -2.5000, -4.0000, 3.0000]
Log of 0 gives ``-inf``::
sage: finance.TimeSeries([1,0,3]).log()[1] -inf """ cdef Py_ssize_t i
def exp(self): """ Return new time series got by applying the exponential map to all the terms in the time series.
OUTPUT:
A new time series.
EXAMPLES::
sage: v = finance.TimeSeries([1..5]); v [1.0000, 2.0000, 3.0000, 4.0000, 5.0000] sage: v.exp() [2.7183, 7.3891, 20.0855, 54.5982, 148.4132] sage: v.exp().log() [1.0000, 2.0000, 3.0000, 4.0000, 5.0000] """ cdef Py_ssize_t i
def abs(self): """ Return new time series got by replacing all entries of ``self`` by their absolute value.
OUTPUT:
A new time series.
EXAMPLES::
sage: v = finance.TimeSeries([1,3.1908439,-4,5.93932]) sage: v [1.0000, 3.1908, -4.0000, 5.9393] sage: v.abs() [1.0000, 3.1908, 4.0000, 5.9393] """ cdef Py_ssize_t i
def diffs(self, Py_ssize_t k = 1): r""" Return the new time series got by taking the differences of successive terms in the time series. So if ``self`` is the time series `X_0, X_1, X_2, \dots`, then this function outputs the series `X_1 - X_0, X_2 - X_1, \dots`. The output series has one less term than the input series. If the optional parameter `k` is given, return `X_k - X_0, X_{2k} - X_k, \dots`.
INPUT:
- ``k`` -- positive integer (default: 1)
OUTPUT:
A new time series.
EXAMPLES::
sage: v = finance.TimeSeries([5,4,1.3,2,8]); v [5.0000, 4.0000, 1.3000, 2.0000, 8.0000] sage: v.diffs() [-1.0000, -2.7000, 0.7000, 6.0000] """ return self.scale_time(k).diffs() cdef Py_ssize_t i
def scale_time(self, Py_ssize_t k): r""" Return the new time series at scale ``k``. If the input time series is `X_0, X_1, X_2, \dots`, then this function returns the shorter time series `X_0, X_k, X_{2k}, \dots`.
INPUT:
- ``k`` -- a positive integer.
OUTPUT:
A new time series.
EXAMPLES::
sage: v = finance.TimeSeries([5,4,1.3,2,8,10,3,-5]); v [5.0000, 4.0000, 1.3000, 2.0000, 8.0000, 10.0000, 3.0000, -5.0000] sage: v.scale_time(1) [5.0000, 4.0000, 1.3000, 2.0000, 8.0000, 10.0000, 3.0000, -5.0000] sage: v.scale_time(2) [5.0000, 1.3000, 8.0000, 3.0000] sage: v.scale_time(3) [5.0000, 2.0000] sage: v.scale_time(10) []
A series of odd length::
sage: v = finance.TimeSeries([1..5]); v [1.0000, 2.0000, 3.0000, 4.0000, 5.0000] sage: v.scale_time(2) [1.0000, 3.0000, 5.0000]
TESTS::
sage: v.scale_time(0) Traceback (most recent call last): ... ValueError: k must be positive sage: v.scale_time(-1) Traceback (most recent call last): ... ValueError: k must be positive """
cdef Py_ssize_t i, n
cpdef rescale(self, double s): """ Change ``self`` by multiplying every value in the series by ``s``.
INPUT:
- ``s`` -- a float.
EXAMPLES::
sage: v = finance.TimeSeries([5,4,1.3,2,8,10,3,-5]); v [5.0000, 4.0000, 1.3000, 2.0000, 8.0000, 10.0000, 3.0000, -5.0000] sage: v.rescale(0.5) sage: v [2.5000, 2.0000, 0.6500, 1.0000, 4.0000, 5.0000, 1.5000, -2.5000] """
def scale(self, double s): """ Return new time series obtained by multiplying every value in the series by ``s``.
INPUT:
- ``s`` -- a float.
OUTPUT:
A new time series with all values multiplied by ``s``.
EXAMPLES::
sage: v = finance.TimeSeries([5,4,1.3,2,8,10,3,-5]); v [5.0000, 4.0000, 1.3000, 2.0000, 8.0000, 10.0000, 3.0000, -5.0000] sage: v.scale(0.5) [2.5000, 2.0000, 0.6500, 1.0000, 4.0000, 5.0000, 1.5000, -2.5000] """
def add_scalar(self, double s): """ Return new time series obtained by adding a scalar to every value in the series.
.. NOTE::
To add componentwise, use the ``add_entries`` method.
INPUT:
- ``s`` -- a float.
OUTPUT:
A new time series with ``s`` added to all values.
EXAMPLES::
sage: v = finance.TimeSeries([5,4,1.3,2,8,10,3,-5]); v [5.0000, 4.0000, 1.3000, 2.0000, 8.0000, 10.0000, 3.0000, -5.0000] sage: v.add_scalar(0.5) [5.5000, 4.5000, 1.8000, 2.5000, 8.5000, 10.5000, 3.5000, -4.5000] """
def add_entries(self, t): """ Add corresponding entries of ``self`` and ``t`` together, extending either ``self`` or ``t`` by 0's if they do not have the same length.
.. NOTE::
To add a single number to the entries of a time series, use the ``add_scalar`` method.
INPUT:
- ``t`` -- a time series.
OUTPUT:
A time series with length the maxima of the lengths of ``self`` and ``t``.
EXAMPLES::
sage: v = finance.TimeSeries([1,2,-5]); v [1.0000, 2.0000, -5.0000] sage: v.add_entries([3,4]) [4.0000, 6.0000, -5.0000] sage: v.add_entries(v) [2.0000, 4.0000, -10.0000] sage: v.add_entries([3,4,7,18.5]) [4.0000, 6.0000, 2.0000, 18.5000] """ cdef Py_ssize_t i, n else:
# Iterate up to the length of the shorter one
# Copy over the rest of the values longer._values + shorter._length, sizeof(double)*(v._length - shorter._length))
def show(self, *args, **kwds): """ Calls plot and passes all arguments onto the plot function. This is thus just an alias for plot.
EXAMPLES:
Draw a plot of a time series::
sage: finance.TimeSeries([1..10]).show() Graphics object consisting of 1 graphics primitive """
def plot(self, Py_ssize_t plot_points=1000, points=False, **kwds): r""" Return a plot of this time series as a line or points through `(i,T(i))`, where `i` ranges over nonnegative integers up to the length of ``self``.
INPUT:
- ``plot_points`` -- (default: 1000) 0 or positive integer. Only plot the given number of equally spaced points in the time series. If 0, plot all points.
- ``points`` -- bool (default: ``False``). If ``True``, return just the points of the time series.
- ``**kwds`` -- passed to the line or point command.
EXAMPLES::
sage: v = finance.TimeSeries([5,4,1.3,2,8,10,3,-5]); v [5.0000, 4.0000, 1.3000, 2.0000, 8.0000, 10.0000, 3.0000, -5.0000] sage: v.plot() Graphics object consisting of 1 graphics primitive sage: v.plot(points=True) Graphics object consisting of 1 graphics primitive sage: v.plot() + v.plot(points=True, rgbcolor='red') Graphics object consisting of 2 graphics primitives sage: v.plot() + v.plot(points=True, rgbcolor='red', pointsize=50) Graphics object consisting of 2 graphics primitives """ cdef Py_ssize_t s
s = self._length/plot_points if plot_points <= 0: raise ValueError("plot_points must be a positive integer") v = self.scale_time(s).list()[:plot_points] else: else:
def simple_moving_average(self, Py_ssize_t k): """ Return the moving average time series over the last ``k`` time units. Assumes the input time series was constant with its starting value for negative time. The t-th step of the output is the sum of the previous ``k - 1`` steps of ``self`` and the ``k``-th step divided by ``k``. Thus ``k`` values are averaged at each point.
INPUT:
- ``k`` -- positive integer.
OUTPUT:
A time series with the same number of steps as ``self``.
EXAMPLES::
sage: v = finance.TimeSeries([1,1,1,2,3]); v [1.0000, 1.0000, 1.0000, 2.0000, 3.0000] sage: v.simple_moving_average(0) [1.0000, 1.0000, 1.0000, 2.0000, 3.0000] sage: v.simple_moving_average(1) [1.0000, 1.0000, 1.0000, 2.0000, 3.0000] sage: v.simple_moving_average(2) [1.0000, 1.0000, 1.0000, 1.5000, 2.5000] sage: v.simple_moving_average(3) [1.0000, 1.0000, 1.0000, 1.3333, 2.0000] """ raise ValueError("k must be positive") cdef Py_ssize_t i return t
else: # I have a slight concern about accumulated rounding error given how # this algorithm adds and subtracts.
def exponential_moving_average(self, double alpha): """ Return the exponential moving average time series. Assumes the input time series was constant with its starting value for negative time. The t-th step of the output is the sum of the previous k-1 steps of ``self`` and the k-th step divided by k.
The 0-th term is formally undefined, so we define it to be 0, and we define the first term to be ``self[0]``.
INPUT:
- ``alpha`` -- float; a smoothing factor with ``0 <= alpha <= 1``.
OUTPUT:
A time series with the same number of steps as ``self``.
EXAMPLES::
sage: v = finance.TimeSeries([1,1,1,2,3]); v [1.0000, 1.0000, 1.0000, 2.0000, 3.0000] sage: v.exponential_moving_average(0) [0.0000, 1.0000, 1.0000, 1.0000, 1.0000] sage: v.exponential_moving_average(1) [0.0000, 1.0000, 1.0000, 1.0000, 2.0000] sage: v.exponential_moving_average(0.5) [0.0000, 1.0000, 1.0000, 1.0000, 1.5000]
Some more examples::
sage: v = finance.TimeSeries([1,2,3,4,5]) sage: v.exponential_moving_average(1) [0.0000, 1.0000, 2.0000, 3.0000, 4.0000] sage: v.exponential_moving_average(0) [0.0000, 1.0000, 1.0000, 1.0000, 1.0000] """ raise ValueError("alpha must be between 0 and 1") cdef Py_ssize_t i return t return t
def sums(self, double s=0): """ Return the new time series got by taking the running partial sums of the terms of this time series.
INPUT:
- ``s`` -- starting value for partial sums.
OUTPUT:
A time series.
EXAMPLES::
sage: v = finance.TimeSeries([1,1,1,2,3]); v [1.0000, 1.0000, 1.0000, 2.0000, 3.0000] sage: v.sums() [1.0000, 2.0000, 3.0000, 5.0000, 8.0000] """ cdef Py_ssize_t i
cpdef double sum(self): """ Return the sum of all the entries of ``self``. If ``self`` has length 0, returns 0.
OUTPUT:
A double.
EXAMPLES::
sage: v = finance.TimeSeries([1,1,1,2,3]); v [1.0000, 1.0000, 1.0000, 2.0000, 3.0000] sage: v.sum() 8.0 """ cdef Py_ssize_t i
def prod(self): """ Return the product of all the entries of ``self``. If ``self`` has length 0, returns 1.
OUTPUT:
A 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 Py_ssize_t i
def mean(self): """ Return the mean (average) of the elements of ``self``.
OUTPUT:
A double.
EXAMPLES::
sage: v = finance.TimeSeries([1,1,1,2,3]); v [1.0000, 1.0000, 1.0000, 2.0000, 3.0000] sage: v.mean() 1.6 """
def pow(self, double k): """ Return a new time series with every elements of ``self`` raised to the k-th power.
INPUT:
- ``k`` -- a float.
OUTPUT:
A time series.
EXAMPLES::
sage: v = finance.TimeSeries([1,1,1,2,3]); v [1.0000, 1.0000, 1.0000, 2.0000, 3.0000] sage: v.pow(2) [1.0000, 1.0000, 1.0000, 4.0000, 9.0000] """ cdef Py_ssize_t i
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:
A 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.6 sage: v.moment(2) 3.2 """ raise ValueError("k must be positive") cdef Py_ssize_t i
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:
A double.
EXAMPLES::
sage: v = finance.TimeSeries([1,2,3]) sage: v.central_moment(2) 0.6666666666666666
Note that the central moment is different from the moment here, since the mean is not `0`::
sage: v.moment(2) # doesn't subtract off mean 4.666666666666667
We compute the central moment directly::
sage: mu = v.mean(); mu 2.0 sage: ((1-mu)^2 + (2-mu)^2 + (3-mu)^2) / 3 0.6666666666666666 """ return float(0) # We could make this slightly faster by doing the add scalar # and moment calculation together. But that would be nasty.
def covariance(self, TimeSeries other): r""" Return the covariance of the time series ``self`` and ``other``.
INPUT:
- ``self``, ``other`` -- time series.
Whichever time series has more terms is truncated.
EXAMPLES::
sage: v = finance.TimeSeries([1,-2,3]); w = finance.TimeSeries([4,5,-10]) sage: v.covariance(w) -11.777777777777779 """ cdef Py_ssize_t i # NOTE: There is also a formula # self._values[i]*other._values[i] - mu*mu2 # but that seems less numerically stable (?), and when tested # was not noticeably faster.
def autocovariance(self, Py_ssize_t k=0): r""" Return the k-th autocovariance function `\gamma(k)` of ``self``. This is the covariance of ``self`` with ``self`` shifted by `k`, i.e.,
.. MATH::
\left. \left( \sum_{t=0}^{n-k-1} (x_t - \mu)(x_{t + k} - \mu) \right) \right/ n,
where `n` is the length of ``self``.
Note the denominator of `n`, which gives a "better" sample estimator.
INPUT:
- ``k`` -- a nonnegative integer (default: 0)
OUTPUT:
A float.
EXAMPLES::
sage: v = finance.TimeSeries([13,8,15,4,4,12,11,7,14,12]) sage: v.autocovariance(0) 14.4 sage: mu = v.mean(); sum([(a-mu)^2 for a in v])/len(v) 14.4 sage: v.autocovariance(1) -2.7 sage: mu = v.mean(); sum([(v[i]-mu)*(v[i+1]-mu) for i in range(len(v)-1)])/len(v) -2.7 sage: v.autocovariance(1) -2.7
We illustrate with a random sample that an independently and identically distributed distribution with zero mean and variance `\sigma^2` has autocovariance function `\gamma(h)` with `\gamma(0) = \sigma^2` and `\gamma(h) = 0` for `h\neq 0`. ::
sage: set_random_seed(0) sage: v = finance.TimeSeries(10^6) sage: v.randomize('normal', 0, 5) [3.3835, -2.0055, 1.7882, -2.9319, -4.6827 ... -5.1868, 9.2613, 0.9274, -6.2282, -8.7652] sage: v.autocovariance(0) 24.95410689... sage: v.autocovariance(1) -0.00508390... sage: v.autocovariance(2) 0.022056325... sage: v.autocovariance(3) -0.01902000... """ cdef Py_ssize_t i
def correlation(self, TimeSeries other): """ Return the correlation of ``self`` and ``other``, which is the covariance of ``self`` and ``other`` divided by the product of their standard deviation.
INPUT:
- ``self``, ``other`` -- time series.
Whichever time series has more terms is truncated.
EXAMPLES::
sage: v = finance.TimeSeries([1,-2,3]); w = finance.TimeSeries([4,5,-10]) sage: v.correlation(w) -0.558041609... sage: v.covariance(w)/(v.standard_deviation() * w.standard_deviation()) -0.558041609... """
def autocorrelation(self, Py_ssize_t k=1): r""" Return the k-th sample autocorrelation of this time series `x_i`.
Let `\mu` be the sample mean. Then the sample autocorrelation function is
.. MATH::
\frac{\sum_{t=0}^{n-k-1} (x_t - \mu)(x_{t+k} - \mu) } {\sum_{t=0}^{n-1} (x_t - \mu)^2}.
Note that the variance must be nonzero or you will get a ``ZeroDivisionError``.
INPUT:
- ``k`` -- a nonnegative integer (default: 1)
OUTPUT:
A time series.
EXAMPLES::
sage: v = finance.TimeSeries([13,8,15,4,4,12,11,7,14,12]) sage: v.autocorrelation() -0.1875 sage: v.autocorrelation(1) -0.1875 sage: v.autocorrelation(0) 1.0 sage: v.autocorrelation(2) -0.20138888888888887 sage: v.autocorrelation(3) 0.18055555555555555
sage: finance.TimeSeries([1..1000]).autocorrelation() 0.997 """
def variance(self, bias=False): """ Return the variance of the elements of ``self``, which is the mean of the squares of the differences from the mean.
INPUT:
- ``bias`` -- bool (default: ``False``); if ``False``, divide by ``self.length() - 1`` instead of ``self.length()`` to give a less biased estimator for the variance.
OUTPUT:
A double.
EXAMPLES::
sage: v = finance.TimeSeries([1,1,1,2,3]); v [1.0000, 1.0000, 1.0000, 2.0000, 3.0000] sage: v.variance() 0.8 sage: v.variance(bias=True) 0.64
TESTS::
sage: finance.TimeSeries([1]).variance() 0.0 sage: finance.TimeSeries([]).variance() 0.0 """ cdef double a cdef Py_ssize_t i else:
def standard_deviation(self, bias=False): """ Return the standard deviation of the entries of ``self``.
INPUT:
- ``bias`` -- bool (default: ``False``); if ``False``, divide by ``self.length() - 1`` instead of ``self.length()`` to give a less biased estimator for the variance.
OUTPUT:
A double.
EXAMPLES::
sage: v = finance.TimeSeries([1,1,1,2,3]); v [1.0000, 1.0000, 1.0000, 2.0000, 3.0000] sage: v.standard_deviation() 0.8944271909... sage: v.standard_deviation(bias=True) 0.8
TESTS::
sage: finance.TimeSeries([1]).standard_deviation() 0.0 sage: finance.TimeSeries([]).standard_deviation() 0.0 """
def range_statistic(self, b=None): r""" Return the rescaled range statistic `R/S` of ``self``, which is defined as follows (see Hurst 1951). If the optional parameter ``b`` is given, return the average of `R/S` range statistics of disjoint blocks of size ``b``.
Let `\sigma` be the standard deviation of the sequence of differences of ``self``, and let `Y_k` be the k-th term of ``self``. Let `n` be the number of terms of ``self``, and set `Z_k = Y_k - ((k+1)/n) \cdot Y_n`. Then
.. MATH::
R/S = \big( \max(Z_k) - \min(Z_k) \big) / \sigma
where the max and min are over all `Z_k`. Basically replacing `Y_k` by `Z_k` allows us to measure the difference from the line from the origin to `(n,Y_n)`.
INPUT:
- ``self`` -- a time series (*not* the series of differences).
- ``b`` -- integer (default: ``None``); if given instead divide the input time series up into ``j = floor(n/b)`` disjoint blocks, compute the `R/S` statistic for each block, and return the average of those `R/S` statistics.
OUTPUT:
A float.
EXAMPLES:
Notice that if we make a Brownian motion random walk, there is no difference if we change the standard deviation. ::
sage: set_random_seed(0); finance.TimeSeries(10^6).randomize('normal').sums().range_statistic() 1897.8392602... sage: set_random_seed(0); finance.TimeSeries(10^6).randomize('normal',0,100).sums().range_statistic() 1897.8392602... """
def hurst_exponent(self): """ Returns an estimate of the Hurst exponent of this time series. We use the algorithm from pages 61 -- 63 of [Peteres, Fractal Market Analysis (1994); see Google Books].
We define the Hurst exponent of a constant time series to be 1.
EXAMPLES:
The Hurst exponent of Brownian motion is 1/2. We approximate it with some specific samples. Note that the estimator is biased and slightly overestimates. ::
sage: set_random_seed(0) sage: bm = finance.TimeSeries(10^5).randomize('normal').sums(); bm [0.6767, 0.2756, 0.6332, 0.0469, -0.8897 ... 152.2437, 151.5327, 152.7629, 152.9169, 152.9084] sage: bm.hurst_exponent() 0.527450972...
We compute the Hurst exponent of a simulated fractional Brownian motion with Hurst parameter 0.7. This function estimates the Hurst exponent as 0.706511951... ::
sage: set_random_seed(0) sage: fbm = finance.fractional_brownian_motion_simulation(0.7,0.1,10^5,1)[0] sage: fbm.hurst_exponent() 0.706511951...
Another example with small Hurst exponent (notice the overestimation).
::
sage: fbm = finance.fractional_brownian_motion_simulation(0.2,0.1,10^5,1)[0] sage: fbm.hurst_exponent() 0.278997441...
We compute the mean Hurst exponent of 100 simulated multifractal cascade random walks::
sage: set_random_seed(0) sage: y = finance.multifractal_cascade_random_walk_simulation(3700,0.02,0.01,0.01,1000,100) sage: finance.TimeSeries([z.hurst_exponent() for z in y]).mean() 0.57984822577934...
We compute the mean Hurst exponent of 100 simulated Markov switching multifractal time series. The Hurst exponent is quite small. ::
sage: set_random_seed(0) sage: msm = finance.MarkovSwitchingMultifractal(8,1.4,0.5,0.95,3) sage: y = msm.simulations(1000,100) sage: finance.TimeSeries([z.hurst_exponent() for z in y]).mean() 0.286102325623705... """ # We use disjoint blocks of size 8, 16, 32, .... k = 4 except ZeroDivisionError: # 0 standard deviation pass return float(1) return v1[0]/v0[0]
def min(self, bint index=False): """ Return the smallest value in this time series. If this series has length 0 we raise a ``ValueError``.
INPUT:
- ``index`` -- bool (default: ``False``); if ``True``, also return index of minimal entry.
OUTPUT:
- float -- smallest value.
- integer -- index of smallest value; only returned if ``index=True``.
EXAMPLES::
sage: v = finance.TimeSeries([1,-4,3,-2.5,-4]) sage: v.min() -4.0 sage: v.min(index=True) (-4.0, 1) """ raise ValueError("min() arg is an empty sequence") cdef Py_ssize_t i, j else:
def max(self, bint index=False): """ Return the largest value in this time series. If this series has length 0 we raise a ``ValueError``.
INPUT:
- ``index`` -- bool (default: ``False``); if ``True``, also return index of maximum entry.
OUTPUT:
- float -- largest value.
- integer -- index of largest value; only returned if ``index=True``.
EXAMPLES::
sage: v = finance.TimeSeries([1,-4,3,-2.5,-4,3]) sage: v.max() 3.0 sage: v.max(index=True) (3.0, 2) """ raise ValueError("max() arg is an empty sequence") else:
def clip_remove(self, min=None, max=None): """ Return new time series obtained from ``self`` by removing all values that are less than or equal to a certain minimum value or greater than or equal to a certain maximum.
INPUT:
- ``min`` -- (default: ``None``) ``None`` or double.
- ``max`` -- (default: ``None``) ``None`` or double.
OUTPUT:
A time series.
EXAMPLES::
sage: v = finance.TimeSeries([1..10]) sage: v.clip_remove(3,7) [3.0000, 4.0000, 5.0000, 6.0000, 7.0000] sage: v.clip_remove(3) [3.0000, 4.0000, 5.0000, 6.0000, 7.0000, 8.0000, 9.0000, 10.0000] sage: v.clip_remove(max=7) [1.0000, 2.0000, 3.0000, 4.0000, 5.0000, 6.0000, 7.0000] """ cdef Py_ssize_t i, j cdef TimeSeries t cdef double mn, mx cdef double x return self.copy() # This code is ugly but I think as fast as possible. # Return everything <= max
# Return everything >= min else: # Return everything between min and max
def histogram(self, Py_ssize_t bins=50, bint normalize=False): """ Return the frequency histogram of the values in this time series divided up into the given number of bins.
INPUT:
- ``bins`` -- a positive integer (default: 50)
- ``normalize`` -- (default: ``False``) whether to normalize so the total area in the bars of the histogram is 1.
OUTPUT:
- counts -- list of counts of numbers of elements in each bin.
- endpoints -- list of 2-tuples (a,b) that give the endpoints of the bins.
EXAMPLES::
sage: v = finance.TimeSeries([5,4,1.3,2,8,10,3,-5]); v [5.0000, 4.0000, 1.3000, 2.0000, 8.0000, 10.0000, 3.0000, -5.0000] sage: v.histogram(3) ([1, 4, 3], [(-5.0, 0.0), (0.0, 5.0), (5.0, 10.0)]) """ raise ValueError("bins must be positive")
cdef Py_ssize_t j
raise ValueError("bins have 0 width")
return [], v
counts = [0]*bins counts[0] = [self._length] return counts, v
cdef Py_ssize_t i raise MemoryError
else:
def plot_histogram(self, bins=50, normalize=True, **kwds): """ Return histogram plot of this time series with given number of bins.
INPUT:
- ``bins`` -- positive integer (default: 50)
- ``normalize`` -- (default: ``True``) whether to normalize so the total area in the bars of the histogram is 1.
OUTPUT:
A histogram plot.
EXAMPLES::
sage: v = finance.TimeSeries([1..50]) sage: v.plot_histogram(bins=10) Graphics object consisting of 10 graphics primitives
::
sage: v.plot_histogram(bins=3,normalize=False,aspect_ratio=1) Graphics object consisting of 3 graphics primitives """
def plot_candlestick(self, int bins=30): """ Return a candlestick plot of this time series with the given number of bins.
A candlestick plot is a style of bar-chart used to view open, high, low, and close stock data. At each bin, the line represents the high / low range. The bar represents the open / close range. The interval is colored blue if the open for that bin is less than the close. If the close is less than the open, then that bin is colored red instead.
INPUT:
- ``bins`` -- positive integer (default: 30), the number of bins or candles.
OUTPUT:
A candlestick plot.
EXAMPLES:
Here we look at the candlestick plot for Brownian motion::
sage: v = finance.TimeSeries(1000).randomize() sage: v.plot_candlestick(bins=20) Graphics object consisting of 40 graphics primitives """
cdef TimeSeries s cdef int i, j, n
def numpy(self, copy=True): """ Return a NumPy version of this time series.
.. NOTE::
If copy is ``False``, return a NumPy 1-D array reference to exactly the same block of memory as this time series. This is very, very fast and makes it easy to quickly use all NumPy/SciPy functionality on ``self``. However, it is dangerous because when this time series goes out of scope and is garbage collected, the corresponding NumPy reference object will point to garbage.
INPUT:
- ``copy`` -- bool (default: ``True``)
OUTPUT:
A numpy 1-D array.
EXAMPLES::
sage: v = finance.TimeSeries([1,-3,4.5,-2]) sage: w = v.numpy(copy=False); w array([ 1. , -3. , 4.5, -2. ]) sage: type(w) <... 'numpy.ndarray'> sage: w.shape (4,)
Notice that changing ``w`` also changes ``v`` too! ::
sage: w[0] = 20 sage: w array([ 20. , -3. , 4.5, -2. ]) sage: v [20.0000, -3.0000, 4.5000, -2.0000]
If you want a separate copy do not give the ``copy=False`` option. ::
sage: z = v.numpy(); z array([ 20. , -3. , 4.5, -2. ]) sage: z[0] = -10 sage: v [20.0000, -3.0000, 4.5000, -2.0000] """ cdef cnumpy.npy_intp dims[1] else: # Py_INCREF(self) # n.base = self
def randomize(self, distribution='uniform', loc=0, scale=1, **kwds): r""" Randomize the entries in this time series, and return a reference to ``self``. Thus this function both changes ``self`` in place, and returns a copy of ``self``, for convenience.
INPUT:
- ``distribution`` -- (default: ``"uniform"``); supported values are:
- ``'uniform'`` -- from ``loc`` to ``loc + scale``
- ``'normal'`` -- mean ``loc`` and standard deviation ``scale``
- ``'semicircle'`` -- with center at ``loc`` (``scale`` is ignored)
- ``'lognormal'`` -- mean ``loc`` and standard deviation ``scale``
- ``loc`` -- float (default: 0)
- ``scale`` -- float (default: 1)
.. NOTE::
All random numbers are generated using algorithms that build on the high quality GMP random number function gmp_urandomb_ui. Thus this function fully respects the Sage ``set_random_state`` command. It's not quite as fast as the C library random number generator, but is of course much better quality, and is platform independent.
EXAMPLES:
We generate 5 uniform random numbers in the interval [0,1]::
sage: set_random_seed(0) sage: finance.TimeSeries(5).randomize() [0.8685, 0.2816, 0.0229, 0.1456, 0.7314]
We generate 5 uniform random numbers from 5 to `5+2=7`::
sage: set_random_seed(0) sage: finance.TimeSeries(10).randomize('uniform', 5, 2) [6.7369, 5.5632, 5.0459, 5.2911, 6.4628, 5.2412, 5.2010, 5.2761, 5.5813, 5.5439]
We generate 5 normal random values with mean 0 and variance 1. ::
sage: set_random_seed(0) sage: finance.TimeSeries(5).randomize('normal') [0.6767, -0.4011, 0.3576, -0.5864, -0.9365]
We generate 10 normal random values with mean 5 and variance 2. ::
sage: set_random_seed(0) sage: finance.TimeSeries(10).randomize('normal', 5, 2) [6.3534, 4.1978, 5.7153, 3.8273, 3.1269, 2.9598, 3.7484, 6.7472, 3.8986, 4.6271]
We generate 5 values using the semicircle distribution. ::
sage: set_random_seed(0) sage: finance.TimeSeries(5).randomize('semicircle') [0.7369, -0.9541, 0.4628, -0.7990, -0.4187]
We generate 1 million normal random values and create a frequency histogram. ::
sage: set_random_seed(0) sage: a = finance.TimeSeries(10^6).randomize('normal') sage: a.histogram(10)[0] [36, 1148, 19604, 130699, 340054, 347870, 137953, 21290, 1311, 35]
We take the above values, and compute the proportion that lie within 1, 2, 3, 5, and 6 standard deviations of the mean (0)::
sage: s = a.standard_deviation() sage: len(a.clip_remove(-s,s))/float(len(a)) 0.683094 sage: len(a.clip_remove(-2*s,2*s))/float(len(a)) 0.954559 sage: len(a.clip_remove(-3*s,3*s))/float(len(a)) 0.997228 sage: len(a.clip_remove(-5*s,5*s))/float(len(a)) 0.999998
There were no "six sigma events"::
sage: len(a.clip_remove(-6*s,6*s))/float(len(a)) 1.0 """ else: raise NotImplementedError
def _randomize_uniform(self, double left, double right): """ Generates a uniform random distribution of doubles between ``left`` and ``right`` and stores values in place.
INPUT:
- ``left`` -- left bound on random distribution.
- ``right`` -- right bound on random distribution.
EXAMPLES:
We generate 5 values distributed with respect to the uniform distribution over the interval [0,1]::
sage: v = finance.TimeSeries(5) sage: set_random_seed(0) sage: v.randomize('uniform') [0.8685, 0.2816, 0.0229, 0.1456, 0.7314]
We now test that the mean is indeed 0.5::
sage: v = finance.TimeSeries(10^6) sage: set_random_seed(0) sage: v.randomize('uniform').mean() 0.50069085... """ raise ValueError("left must be less than right")
cdef Py_ssize_t k
def _randomize_normal(self, double m, double s): """ Generates a normal random distribution of doubles with mean ``m`` and standard deviation ``s`` and stores values in place. Uses the Box-Muller algorithm.
INPUT:
- ``m`` -- mean
- ``s`` -- standard deviation
EXAMPLES:
We generate 5 values distributed with respect to the normal distribution with mean 0 and standard deviation 1::
sage: set_random_seed(0) sage: v = finance.TimeSeries(5) sage: v.randomize('normal') [0.6767, -0.4011, 0.3576, -0.5864, -0.9365]
We now test that the mean is indeed 0::
sage: set_random_seed(0) sage: v = finance.TimeSeries(10^6) sage: v.randomize('normal').mean() 6.2705472723...
The same test with mean equal to 2 and standard deviation equal to 5::
sage: set_random_seed(0) sage: v = finance.TimeSeries(10^6) sage: v.randomize('normal', 2, 5).mean() 2.0003135273... """ # Ported from http://users.tkk.fi/~nbeijar/soft/terrain/source_o2/boxmuller.c # This the box muller algorithm.
cdef double x1, x2, w, y1, y2 cdef Py_ssize_t k
def _randomize_semicircle(self, double center): """ Generates a semicircle random distribution of doubles about center and stores values in place. Uses the acceptance-rejection method.
INPUT:
- ``center`` -- the center of the semicircle distribution.
EXAMPLES:
We generate 5 values distributed with respect to the semicircle distribution located at center::
sage: v = finance.TimeSeries(5) sage: set_random_seed(0) sage: v.randomize('semicircle') [0.7369, -0.9541, 0.4628, -0.7990, -0.4187]
We now test that the mean is indeed the center::
sage: v = finance.TimeSeries(10^6) sage: set_random_seed(0) sage: v.randomize('semicircle').mean() 0.0007207497...
The same test with center equal to 2::
sage: v = finance.TimeSeries(10^6) sage: set_random_seed(0) sage: v.randomize('semicircle', 2).mean() 2.0007207497... """ cdef Py_ssize_t k
def _randomize_lognormal(self, double m, double s): r""" Generates a log-normal random distribution of doubles with mean ``m`` and standard deviation ``s``. Uses Box-Muller algorithm and the identity: if `Y` is a random variable with normal distribution then `X = \exp(Y)` is a random variable with log-normal distribution.
INPUT:
- ``m`` -- mean
- ``s`` -- standard deviation
EXAMPLES:
We generate 5 values distributed with respect to the lognormal distribution with mean 0 and standard deviation 1::
sage: set_random_seed(0) sage: v = finance.TimeSeries(5) sage: v.randomize('lognormal') [1.9674, 0.6696, 1.4299, 0.5563, 0.3920]
We now test that the mean is indeed `\sqrt{e}`::
sage: set_random_seed(0) sage: v = finance.TimeSeries(10^6) sage: v.randomize('lognormal').mean() 1.647351973... sage: exp(0.5) 1.648721270...
A log-normal distribution can be simply thought of as the logarithm of a normally distributed data set. We test that here by generating 5 values distributed with respect to the normal distribution with mean 0 and standard deviation 1::
sage: set_random_seed(0) sage: w = finance.TimeSeries(5) sage: w.randomize('normal') [0.6767, -0.4011, 0.3576, -0.5864, -0.9365] sage: exp(w) [1.9674, 0.6696, 1.4299, 0.5563, 0.3920] """ # Ported from http://users.tkk.fi/~nbeijar/soft/terrain/source_o2/boxmuller.c # This the box muller algorithm. cdef double x1, x2, w, y1, y2 cdef Py_ssize_t k
def fft(self, bint overwrite=False): r""" Return the real discrete fast Fourier transform of ``self``, as a real time series:
.. MATH::
[y(0),\Re(y(1)),\Im(y(1)),\dots,\Re(y(n/2))] \text{ if $n$ is even}
[y(0),\Re(y(1)),\Im(y(1)),\dots,\Re(y(n/2)),\Im(y(n/2))] \text{ if $n$ is odd}
where
.. MATH::
y(j) = \sum_{k=0}^{n-1} x[k] \cdot \exp(-\sqrt{-1} \cdot jk \cdot 2\pi/n)
for `j = 0,\dots,n-1`. Note that `y(-j) = y(n-j)`.
EXAMPLES::
sage: v = finance.TimeSeries([1..9]); v [1.0000, 2.0000, 3.0000, 4.0000, 5.0000, 6.0000, 7.0000, 8.0000, 9.0000] sage: w = v.fft(); w [45.0000, -4.5000, 12.3636, -4.5000, 5.3629, -4.5000, 2.5981, -4.5000, 0.7935]
We get just the series of real parts of ::
sage: finance.TimeSeries([w[0]]) + w[1:].scale_time(2) [45.0000, -4.5000, -4.5000, -4.5000, -4.5000]
An example with an even number of terms::
sage: v = finance.TimeSeries([1..10]); v [1.0000, 2.0000, 3.0000, 4.0000, 5.0000, 6.0000, 7.0000, 8.0000, 9.0000, 10.0000] sage: w = v.fft(); w [55.0000, -5.0000, 15.3884, -5.0000, 6.8819, -5.0000, 3.6327, -5.0000, 1.6246, -5.0000] sage: v.fft().ifft() [1.0000, 2.0000, 3.0000, 4.0000, 5.0000, 6.0000, 7.0000, 8.0000, 9.0000, 10.0000] """ y = self else:
def ifft(self, bint overwrite=False): r""" Return the real discrete inverse fast Fourier transform of ``self``, which is also a real time series.
This is the inverse of ``fft()``.
The returned real array contains
.. MATH::
[y(0),y(1),\dots,y(n-1)]
where for `n` is even
.. MATH::
y(j) = 1/n \left( \sum_{k=1}^{n/2-1} (x[2k-1]+\sqrt{-1} \cdot x[2k]) \cdot \exp(\sqrt{-1} \cdot jk \cdot 2pi/n) + c.c. + x[0] + (-1)^j x[n-1] \right)
and for `n` is odd
.. MATH::
y(j) = 1/n \left( \sum_{k=1}^{(n-1)/2} (x[2k-1]+\sqrt{-1} \cdot x[2k]) \cdot \exp(\sqrt{-1} \cdot jk \cdot 2pi/n) + c.c. + x[0] \right)
where `c.c.` denotes complex conjugate of preceding expression.
EXAMPLES::
sage: v = finance.TimeSeries([1..10]); v [1.0000, 2.0000, 3.0000, 4.0000, 5.0000, 6.0000, 7.0000, 8.0000, 9.0000, 10.0000] sage: v.ifft() [5.1000, -5.6876, 1.4764, -1.0774, 0.4249, -0.1000, -0.2249, 0.6663, -1.2764, 1.6988] sage: v.ifft().fft() [1.0000, 2.0000, 3.0000, 4.0000, 5.0000, 6.0000, 7.0000, 8.0000, 9.0000, 10.0000] """ else:
cdef new_time_series(Py_ssize_t length): """ Return a new uninitialized time series of the given length. The entries of the time series are garbage.
INPUT:
- ``length`` -- integer
OUTPUT:
A time series.
EXAMPLES:
This uses ``new_time_series`` implicitly::
sage: v = finance.TimeSeries([1,-3,4.5,-2]) sage: v.__copy__() [1.0000, -3.0000, 4.5000, -2.0000] """ raise ValueError("length must be nonnegative")
""" Version 1 unpickle method.
INPUT:
- ``v`` -- a raw char buffer
EXAMPLES::
sage: v = finance.TimeSeries([1,2,3]) sage: s = v.__reduce__()[1][0] sage: type(s) # py2 <type 'str'> sage: type(s) # py3 <type 'bytes'> sage: sage.finance.time_series.unpickle_time_series_v1(s,3) [1.0000, 2.0000, 3.0000] sage: sage.finance.time_series.unpickle_time_series_v1(s+s,6) [1.0000, 2.0000, 3.0000, 1.0000, 2.0000, 3.0000] sage: sage.finance.time_series.unpickle_time_series_v1(b'',0) [] """
r""" Given a sequence of lagged autocovariances of length `M` produce `a_1,\dots,a_p` so that the first `M` autocovariance coefficients of the autoregressive processes `X_t=a_1X_{t_1}+\cdots+a_pX_{t-p}+Z_t` are the same as the input sequence.
The function works by solving the Yule-Walker equations `\Gamma a =\gamma`, where `\gamma=(\gamma(1),\dots,\gamma(M))`, `a=(a_1,\dots,a_M)`, with `\gamma(i)` the autocovariance of lag `i` and `\Gamma_{ij}=\gamma(i-j)`.
EXAMPLES:
In this example we consider the multifractal cascade random walk of length 1000, and use simulations to estimate the expected first few autocovariance parameters for this model, then use them to construct a linear filter that works vastly better than a linear filter constructed from the same data but not using this model. The Monte-Carlo method illustrated below should work for predicting one "time step" into the future for any model that can be simulated. To predict k time steps into the future would require using a similar technique but would require scaling time by k.
We create 100 simulations of a multifractal random walk. This models the logarithms of a stock price sequence. ::
sage: set_random_seed(0) sage: y = finance.multifractal_cascade_random_walk_simulation(3700,0.02,0.01,0.01,1000,100)
For each walk below we replace the walk by the walk but where each step size is replaced by its absolute value -- this is what we expect to be able to predict given the model, which is only a model for predicting volatility. We compute the first 200 autocovariance values for every random walk::
sage: c = [[a.diffs().abs().sums().autocovariance(i) for a in y] for i in range(200)]
We make a time series out of the expected values of the autocovariances::
sage: ac = finance.TimeSeries([finance.TimeSeries(z).mean() for z in c]) sage: ac [3.9962, 3.9842, 3.9722, 3.9601, 3.9481 ... 1.7144, 1.7033, 1.6922, 1.6812, 1.6701]
.. NOTE::
``ac`` looks like a line -- one could best fit it to yield a lot more approximate autocovariances.
We compute the autoregression coefficients matching the above autocovariances::
sage: F = finance.autoregressive_fit(ac); F [0.9982, -0.0002, -0.0002, 0.0003, 0.0001 ... 0.0002, -0.0002, -0.0000, -0.0002, -0.0014]
Note that the sum is close to 1::
sage: sum(F) 0.99593284089454...
Now we make up an 'out of sample' sequence::
sage: y2 = finance.multifractal_cascade_random_walk_simulation(3700,0.02,0.01,0.01,1000,1)[0].diffs().abs().sums() sage: y2 [0.0013, 0.0059, 0.0066, 0.0068, 0.0184 ... 6.8004, 6.8009, 6.8063, 6.8090, 6.8339]
And we forecast the very last value using our linear filter; the forecast is close::
sage: y2[:-1].autoregressive_forecast(F) 6.7836741372407...
In fact it is closer than we would get by forecasting using a linear filter made from all the autocovariances of our sequence::
sage: y2[:-1].autoregressive_forecast(y2[:-1].autoregressive_fit(len(y2))) 6.770168705668...
We record the last 20 forecasts, always using all correct values up to the one we are forecasting::
sage: s1 = sum([(y2[:-i].autoregressive_forecast(F)-y2[-i])^2 for i in range(1,20)])
We do the same, but using the autocovariances of the sample sequence::
sage: F2 = y2[:-100].autoregressive_fit(len(F)) sage: s2 = sum([(y2[:-i].autoregressive_forecast(F2)-y2[-i])^2 for i in range(1,20)])
Our model gives us something that is 15 percent better in this case::
sage: s2/s1 1.15464636102...
How does it compare overall? To find out we do 100 simulations and for each we compute the percent that our model beats naively using the autocovariances of the sample::
sage: y_out = finance.multifractal_cascade_random_walk_simulation(3700,0.02,0.01,0.01,1000,100) sage: s1 = []; s2 = [] sage: for v in y_out: ....: s1.append(sum([(v[:-i].autoregressive_forecast(F)-v[-i])^2 for i in range(1,20)])) ....: F2 = v[:-len(F)].autoregressive_fit(len(F)) ....: s2.append(sum([(v[:-i].autoregressive_forecast(F2)-v[-i])^2 for i in range(1,20)]))
We find that overall the model beats naive linear forecasting by 35 percent! ::
sage: s = finance.TimeSeries([s2[i]/s1[i] for i in range(len(s1))]) sage: s.mean() 1.354073591877... """ cdef TimeSeries c cdef Py_ssize_t i
raise ValueError("M must be positive")
else:
# Also create the numpy vector v with entries c(1), ..., c(M).
# 2. Create the autocovariances numpy 2d array A whose i,j entry # is c(|i-j|).
# 3. Solve the equation A * x = v
# 4. The entries of x give the linear filter coefficients. |