-- SPDX-License-Identifier: MPL-2.0

{- | Timestamped finite observations. Times are caller-defined numeric units,
strictly increasing, with explicit missing values. No implicit calendar, time
zone, sorting, duplicate aggregation, or extrapolation occurs. Reference
operations accept at most 16,384 observations; quadratic operations impose
tighter bounds. A regular series is a timestamped series with a positive step.
-}
module HMath.TimeSeries.Core (
    Series,
    series,
    regularSeries,
    observations,
    seriesValues,
    isRegular,
    Alignment (..),
    align,
    Fill (..),
    fillMissing,
    Resampling (..),
    resample,
    Rolling (..),
    rolling,
    exponentialAverage,
    weightedAverage,
    cumulativeReturns,
    logarithmicReturns,
    difference,
    autocorrelation,
    crossCorrelation,
    durbinLevinson,
    yuleWalker,
    partialAutocorrelation,
    finite,
    require,
    checked,
    validate,
    average,
) where

import qualified Data.Map.Strict as M
import Data.Maybe (catMaybes)
import HMath.Error

-- | Opaque, validated timestamps paired with optional observations.
newtype Series = Series [(Double, Maybe Double)] deriving (Series -> Series -> Bool
(Series -> Series -> Bool)
-> (Series -> Series -> Bool) -> Eq Series
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Series -> Series -> Bool
== :: Series -> Series -> Bool
$c/= :: Series -> Series -> Bool
/= :: Series -> Series -> Bool
Eq, Int -> Series -> ShowS
[Series] -> ShowS
Series -> String
(Int -> Series -> ShowS)
-> (Series -> String) -> ([Series] -> ShowS) -> Show Series
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Series -> ShowS
showsPrec :: Int -> Series -> ShowS
$cshow :: Series -> String
show :: Series -> String
$cshowList :: [Series] -> ShowS
showList :: [Series] -> ShowS
Show)

-- | Test whether Double is neither NaN nor infinite.
finite :: Double -> Bool
finite :: Double -> Bool
finite Double
x = Bool -> Bool
not (Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x)

{- | Given predicate then message, return Right () or time-series-prefixed Left
InvalidArgument.
-}
require :: Bool -> String -> Either MathError ()
require :: Bool -> String -> Either MathError ()
require Bool
p String
s = if Bool
p then () -> Either MathError ()
forall a b. b -> Either a b
Right () else MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument (String
"time series: " String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
s))

{- | Return list unchanged when every sample is finite; otherwise return Left
NonFiniteResult.
-}
checked :: [Double] -> Either MathError [Double]
checked :: [Double] -> Either MathError [Double]
checked [Double]
xs = if (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
xs then [Double] -> Either MathError [Double]
forall a b. b -> Either a b
Right [Double]
xs else MathError -> Either MathError [Double]
forall a b. a -> Either a b
Left MathError
NonFiniteResult

{- | Given minimum sample count then list, require length between that minimum and 16384 and
all finite values. Invalid input returns Left InvalidArgument.
-}
validate :: Int -> [Double] -> Either MathError ()
validate :: Int -> [Double] -> Either MathError ()
validate Int
n [Double]
xs =
    Bool -> String -> Either MathError ()
require
        ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
n Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
16384 Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
xs)
        String
"finite sample count outside supported range"

{- | Unchecked arithmetic average using per-term division by list length; empty list returns
zero by empty sum. Checked statistical callers validate nonempty input first.
-}
average :: [Double] -> Double
average :: [Double] -> Double
average [Double]
xs = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs)) [Double]
xs)

-- | Preserve order and missing values; reject duplicate or nonfinite times.
series :: [(Double, Maybe Double)] -> Either MathError Series
series :: [(Double, Maybe Double)] -> Either MathError Series
series [(Double, Maybe Double)]
xs = do
    Bool -> String -> Either MathError ()
require ([(Double, Maybe Double)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(Double, Maybe Double)]
xs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
16384) String
"too many observations"
    Bool -> String -> Either MathError ()
require (((Double, Maybe Double) -> Bool)
-> [(Double, Maybe Double)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Double -> Bool
finite (Double -> Bool)
-> ((Double, Maybe Double) -> Double)
-> (Double, Maybe Double)
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double, Maybe Double) -> Double
forall a b. (a, b) -> a
fst) [(Double, Maybe Double)]
xs Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite ([Maybe Double] -> [Double]
forall a. [Maybe a] -> [a]
catMaybes (((Double, Maybe Double) -> Maybe Double)
-> [(Double, Maybe Double)] -> [Maybe Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double, Maybe Double) -> Maybe Double
forall a b. (a, b) -> b
snd [(Double, Maybe Double)]
xs))) String
"nonfinite observation"
    Bool -> String -> Either MathError ()
require ([Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and ((Double -> Double -> Bool) -> [Double] -> [Double] -> [Bool]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
(<) (((Double, Maybe Double) -> Double)
-> [(Double, Maybe Double)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double, Maybe Double) -> Double
forall a b. (a, b) -> a
fst [(Double, Maybe Double)]
xs) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 (((Double, Maybe Double) -> Double)
-> [(Double, Maybe Double)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double, Maybe Double) -> Double
forall a b. (a, b) -> a
fst [(Double, Maybe Double)]
xs)))) String
"timestamps must strictly increase"
    Series -> Either MathError Series
forall a b. b -> Either a b
Right ([(Double, Maybe Double)] -> Series
Series [(Double, Maybe Double)]
xs)

-- | Sample times start+i*step. Rounded duplicate times are rejected as well.
regularSeries :: Double -> Double -> [Maybe Double] -> Either MathError Series
regularSeries :: Double -> Double -> [Maybe Double] -> Either MathError Series
regularSeries Double
start Double
step [Maybe Double]
xs = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
start Bool -> Bool -> Bool
&& Double -> Bool
finite Double
step Bool -> Bool -> Bool
&& Double
step Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"finite start and positive step required"
    [(Double, Maybe Double)] -> Either MathError Series
series ([Double] -> [Maybe Double] -> [(Double, Maybe Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double
start Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
step | Int
i <- [Int
0 :: Int .. [Maybe Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Maybe Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] [Maybe Double]
xs)

{- | Extract time-ordered (timestamp,Maybe value) observations from a validated Series.
Missing values remain Nothing.
-}
observations :: Series -> [(Double, Maybe Double)]
observations :: Series -> [(Double, Maybe Double)]
observations (Series [(Double, Maybe Double)]
xs) = [(Double, Maybe Double)]
xs

-- | Complete-data conversion never silently drops a missing observation.
seriesValues :: Series -> Either MathError [Double]
seriesValues :: Series -> Either MathError [Double]
seriesValues (Series [(Double, Maybe Double)]
xs) = ((Double, Maybe Double) -> Either MathError Double)
-> [(Double, Maybe Double)] -> Either MathError [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM (Either MathError Double
-> (Double -> Either MathError Double)
-> Maybe Double
-> Either MathError Double
forall b a. b -> (a -> b) -> Maybe a -> b
maybe (MathError -> Either MathError Double
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"time series: missing observation")) Double -> Either MathError Double
forall a b. b -> Either a b
Right (Maybe Double -> Either MathError Double)
-> ((Double, Maybe Double) -> Maybe Double)
-> (Double, Maybe Double)
-> Either MathError Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double, Maybe Double) -> Maybe Double
forall a b. (a, b) -> b
snd) [(Double, Maybe Double)]
xs

-- | Relative spacing tolerance; fewer than three observations are regular.
isRegular :: Double -> Series -> Either MathError Bool
isRegular :: Double -> Series -> Either MathError Bool
isRegular Double
tolerance (Series [(Double, Maybe Double)]
xs) = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
tolerance Bool -> Bool -> Bool
&& Double
tolerance Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0) String
"invalid spacing tolerance"
    let steps :: [Double]
steps = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 (((Double, Maybe Double) -> Double)
-> [(Double, Maybe Double)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double, Maybe Double) -> Double
forall a b. (a, b) -> a
fst [(Double, Maybe Double)]
xs)) (((Double, Maybe Double) -> Double)
-> [(Double, Maybe Double)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double, Maybe Double) -> Double
forall a b. (a, b) -> a
fst [(Double, Maybe Double)]
xs)
    _ <- [Double] -> Either MathError [Double]
checked [Double]
steps
    Right
        (case steps of [] -> Bool
True; Double
d : [Double]
ds -> (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Double
x -> Double -> Double
forall a. Num a => a -> a
abs (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
d) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
tolerance Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (Double -> Double
forall a. Num a => a -> a
abs Double
d) (Double -> Double
forall a. Num a => a -> a
abs Double
x)) [Double]
ds)

{- | Time-series timestamp join policy: Inner intersection or Outer union,with missing
values explicit.
-}
data Alignment = Inner | Outer | LeftAligned deriving (Alignment -> Alignment -> Bool
(Alignment -> Alignment -> Bool)
-> (Alignment -> Alignment -> Bool) -> Eq Alignment
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Alignment -> Alignment -> Bool
== :: Alignment -> Alignment -> Bool
$c/= :: Alignment -> Alignment -> Bool
/= :: Alignment -> Alignment -> Bool
Eq, Int -> Alignment -> ShowS
[Alignment] -> ShowS
Alignment -> String
(Int -> Alignment -> ShowS)
-> (Alignment -> String)
-> ([Alignment] -> ShowS)
-> Show Alignment
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Alignment -> ShowS
showsPrec :: Int -> Alignment -> ShowS
$cshow :: Alignment -> String
show :: Alignment -> String
$cshowList :: [Alignment] -> ShowS
showList :: [Alignment] -> ShowS
Show)

-- | Exact timestamp equality. Missing alignment positions are explicit Nothing.
align :: Alignment -> Series -> Series -> Either MathError (Series, Series)
align :: Alignment -> Series -> Series -> Either MathError (Series, Series)
align Alignment
mode (Series [(Double, Maybe Double)]
a) (Series [(Double, Maybe Double)]
b) = do
    let am :: Map Double (Maybe Double)
am = [(Double, Maybe Double)] -> Map Double (Maybe Double)
forall k a. Ord k => [(k, a)] -> Map k a
M.fromList [(Double, Maybe Double)]
a
        bm :: Map Double (Maybe Double)
bm = [(Double, Maybe Double)] -> Map Double (Maybe Double)
forall k a. Ord k => [(k, a)] -> Map k a
M.fromList [(Double, Maybe Double)]
b
        times :: [Double]
times = case Alignment
mode of
            Alignment
Inner -> Map Double (Maybe Double) -> [Double]
forall k a. Map k a -> [k]
M.keys (Map Double (Maybe Double)
-> Map Double (Maybe Double) -> Map Double (Maybe Double)
forall k a b. Ord k => Map k a -> Map k b -> Map k a
M.intersection Map Double (Maybe Double)
am Map Double (Maybe Double)
bm)
            Alignment
Outer -> Map Double (Maybe Double) -> [Double]
forall k a. Map k a -> [k]
M.keys (Map Double (Maybe Double)
-> Map Double (Maybe Double) -> Map Double (Maybe Double)
forall k a. Ord k => Map k a -> Map k a -> Map k a
M.union Map Double (Maybe Double)
am Map Double (Maybe Double)
bm)
            Alignment
LeftAligned -> Map Double (Maybe Double) -> [Double]
forall k a. Map k a -> [k]
M.keys Map Double (Maybe Double)
am
    x <- [(Double, Maybe Double)] -> Either MathError Series
series [(Double
t, Maybe Double -> Double -> Map Double (Maybe Double) -> Maybe Double
forall k a. Ord k => a -> k -> Map k a -> a
M.findWithDefault Maybe Double
forall a. Maybe a
Nothing Double
t Map Double (Maybe Double)
am) | Double
t <- [Double]
times]
    y <- series [(t, M.findWithDefault Nothing t bm) | t <- times]
    Right (x, y)

{- | Missing-value strategy selecting forward fill,backward fill,or linear interpolation
under the checked boundary rules.
-}
data Fill = ForwardFill | BackwardFill | LinearFill deriving (Fill -> Fill -> Bool
(Fill -> Fill -> Bool) -> (Fill -> Fill -> Bool) -> Eq Fill
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Fill -> Fill -> Bool
== :: Fill -> Fill -> Bool
$c/= :: Fill -> Fill -> Bool
/= :: Fill -> Fill -> Bool
Eq, Int -> Fill -> ShowS
[Fill] -> ShowS
Fill -> String
(Int -> Fill -> ShowS)
-> (Fill -> String) -> ([Fill] -> ShowS) -> Show Fill
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Fill -> ShowS
showsPrec :: Int -> Fill -> ShowS
$cshow :: Fill -> String
show :: Fill -> String
$cshowList :: [Fill] -> ShowS
showList :: [Fill] -> ShowS
Show)

{- | Directional filling leaves the unavailable leading/trailing end missing.
Linear filling interpolates between finite anchors and never extrapolates.
-}
fillMissing :: Fill -> Series -> Either MathError Series
fillMissing :: Fill -> Series -> Either MathError Series
fillMissing Fill
mode (Series [(Double, Maybe Double)]
xs) =
    [(Double, Maybe Double)] -> Either MathError Series
series
        ( case Fill
mode of
            Fill
ForwardFill -> [(Double, Maybe Double)] -> [(Double, Maybe Double)]
forall {a} {a}. [(a, Maybe a)] -> [(a, Maybe a)]
directional [(Double, Maybe Double)]
xs
            Fill
BackwardFill -> [(Double, Maybe Double)] -> [(Double, Maybe Double)]
forall a. [a] -> [a]
reverse ([(Double, Maybe Double)] -> [(Double, Maybe Double)]
forall {a} {a}. [(a, Maybe a)] -> [(a, Maybe a)]
directional ([(Double, Maybe Double)] -> [(Double, Maybe Double)]
forall a. [a] -> [a]
reverse [(Double, Maybe Double)]
xs))
            Fill
LinearFill -> ((Double, Maybe Double) -> (Double, Maybe Double))
-> [(Double, Maybe Double)] -> [(Double, Maybe Double)]
forall a b. (a -> b) -> [a] -> [b]
map (Double, Maybe Double) -> (Double, Maybe Double)
interpolate [(Double, Maybe Double)]
xs
        )
  where
    directional :: [(a, Maybe a)] -> [(a, Maybe a)]
directional = Maybe a -> [(a, Maybe a)] -> [(a, Maybe a)]
forall {a} {a}. Maybe a -> [(a, Maybe a)] -> [(a, Maybe a)]
go Maybe a
forall a. Maybe a
Nothing
    go :: Maybe a -> [(a, Maybe a)] -> [(a, Maybe a)]
go Maybe a
_ [] = []
    go Maybe a
old ((a
t, Maybe a
v) : [(a, Maybe a)]
rest) = let next :: Maybe a
next = Maybe a -> (a -> Maybe a) -> Maybe a -> Maybe a
forall b a. b -> (a -> b) -> Maybe a -> b
maybe Maybe a
old a -> Maybe a
forall a. a -> Maybe a
Just Maybe a
v in (a
t, Maybe a
next) (a, Maybe a) -> [(a, Maybe a)] -> [(a, Maybe a)]
forall a. a -> [a] -> [a]
: Maybe a -> [(a, Maybe a)] -> [(a, Maybe a)]
go Maybe a
next [(a, Maybe a)]
rest
    anchors :: Map Double Double
anchors = [(Double, Double)] -> Map Double Double
forall k a. Ord k => [(k, a)] -> Map k a
M.fromList [(Double
t, Double
v) | (Double
t, Just Double
v) <- [(Double, Maybe Double)]
xs]
    interpolate :: (Double, Maybe Double) -> (Double, Maybe Double)
interpolate pair :: (Double, Maybe Double)
pair@(Double
_, Just Double
_) = (Double, Maybe Double)
pair
    interpolate pair :: (Double, Maybe Double)
pair@(Double
t, Maybe Double
Nothing) = case (Double -> Map Double Double -> Maybe (Double, Double)
forall k v. Ord k => k -> Map k v -> Maybe (k, v)
M.lookupLT Double
t Map Double Double
anchors, Double -> Map Double Double -> Maybe (Double, Double)
forall k v. Ord k => k -> Map k v -> Maybe (k, v)
M.lookupGT Double
t Map Double Double
anchors) of
        (Just (Double
a, Double
x), Just (Double
b, Double
y)) -> (Double
t, Double -> Maybe Double
forall a. a -> Maybe a
Just ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y)) where w :: Double
w = Double -> Double -> Double -> Double
fraction Double
a Double
b Double
t
        (Maybe (Double, Double), Maybe (Double, Double))
_ -> (Double, Maybe Double)
pair

-- | Series resampling policy selecting bin means,last-bin values,or linear point samples.
data Resampling = LinearSamples | BinMean | BinLast deriving (Resampling -> Resampling -> Bool
(Resampling -> Resampling -> Bool)
-> (Resampling -> Resampling -> Bool) -> Eq Resampling
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Resampling -> Resampling -> Bool
== :: Resampling -> Resampling -> Bool
$c/= :: Resampling -> Resampling -> Bool
/= :: Resampling -> Resampling -> Bool
Eq, Int -> Resampling -> ShowS
[Resampling] -> ShowS
Resampling -> String
(Int -> Resampling -> ShowS)
-> (Resampling -> String)
-> ([Resampling] -> ShowS)
-> Show Resampling
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Resampling -> ShowS
showsPrec :: Int -> Resampling -> ShowS
$cshow :: Resampling -> String
show :: Resampling -> String
$cshowList :: [Resampling] -> ShowS
showList :: [Resampling] -> ShowS
Show)

{- | Linear samples require brackets (or an exact timestamp). Bins are
[start+i*step,start+(i+1)*step); the final right endpoint is excluded.
-}
resample :: Resampling -> Double -> Double -> Int -> Series -> Either MathError Series
resample :: Resampling
-> Double -> Double -> Int -> Series -> Either MathError Series
resample Resampling
method Double
start Double
step Int
count input :: Series
input@(Series [(Double, Maybe Double)]
xs) = do
    Bool -> String -> Either MathError ()
require
        (Int
count Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
count Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
16384 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
start Bool -> Bool -> Bool
&& Double -> Bool
finite Double
step Bool -> Bool -> Bool
&& Double
step Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0)
        String
"invalid resampling grid"
    let targets :: [Double]
targets = [Double
start Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
step | Int
i <- [Int
0 .. Int
count Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
    _ <- [Double] -> Either MathError [Double]
checked ([Double]
targets [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double
start Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
count Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
step | Resampling
method Resampling -> Resampling -> Bool
forall a. Eq a => a -> a -> Bool
/= Resampling
LinearSamples])
    case method of
        Resampling
LinearSamples -> do
            let anchors :: Map Double (Maybe Double)
anchors = [(Double, Maybe Double)] -> Map Double (Maybe Double)
forall k a. Ord k => [(k, a)] -> Map k a
M.fromList [(Double, Maybe Double)]
xs
                lookupValue :: Double -> Maybe Double
lookupValue Double
t = case Double -> Map Double (Maybe Double) -> Maybe (Maybe Double)
forall k a. Ord k => k -> Map k a -> Maybe a
M.lookup Double
t Map Double (Maybe Double)
anchors of
                    Just Maybe Double
v -> Maybe Double
v
                    Maybe (Maybe Double)
Nothing -> case (Double -> Map Double (Maybe Double) -> Maybe (Double, Maybe Double)
forall k v. Ord k => k -> Map k v -> Maybe (k, v)
M.lookupLT Double
t Map Double (Maybe Double)
anchors, Double -> Map Double (Maybe Double) -> Maybe (Double, Maybe Double)
forall k v. Ord k => k -> Map k v -> Maybe (k, v)
M.lookupGT Double
t Map Double (Maybe Double)
anchors) of
                        (Just (Double
a, Just Double
x), Just (Double
b, Just Double
y)) -> Double -> Maybe Double
forall a. a -> Maybe a
Just ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y) where w :: Double
w = Double -> Double -> Double -> Double
fraction Double
a Double
b Double
t
                        (Maybe (Double, Maybe Double), Maybe (Double, Maybe Double))
_ -> Maybe Double
forall a. Maybe a
Nothing
            [(Double, Maybe Double)] -> Either MathError Series
series [(Double
t, Double -> Maybe Double
lookupValue Double
t) | Double
t <- [Double]
targets]
        Resampling
_ -> do
            Bool -> String -> Either MathError ()
require
                (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
count Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([(Double, Maybe Double)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Series -> [(Double, Maybe Double)]
observations Series
input)) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
4000000)
                String
"resampling work limit"
            [(Double, Maybe Double)] -> Either MathError Series
series [(Double
t, [Double] -> Maybe Double
aggregate ([Maybe Double] -> [Double]
forall a. [Maybe a] -> [a]
catMaybes [Maybe Double
v | (Double
s, Maybe Double
v) <- [(Double, Maybe Double)]
xs, Double
s Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
t, Double
s Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
step])) | Double
t <- [Double]
targets]
  where
    aggregate :: [Double] -> Maybe Double
aggregate [] = Maybe Double
forall a. Maybe a
Nothing
    aggregate [Double]
vs = Double -> Maybe Double
forall a. a -> Maybe a
Just (if Resampling
method Resampling -> Resampling -> Bool
forall a. Eq a => a -> a -> Bool
== Resampling
BinLast then [Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
vs else [Double] -> Double
average [Double]
vs)

{- | Rolling-window statistic selector for mean,variance,standard deviation,minimum,or
maximum.
-}
data Rolling = MovingMean | MovingVariance | MovingSD | MovingMinimum | MovingMaximum
    deriving (Rolling -> Rolling -> Bool
(Rolling -> Rolling -> Bool)
-> (Rolling -> Rolling -> Bool) -> Eq Rolling
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Rolling -> Rolling -> Bool
== :: Rolling -> Rolling -> Bool
$c/= :: Rolling -> Rolling -> Bool
/= :: Rolling -> Rolling -> Bool
Eq, Int -> Rolling -> ShowS
[Rolling] -> ShowS
Rolling -> String
(Int -> Rolling -> ShowS)
-> (Rolling -> String) -> ([Rolling] -> ShowS) -> Show Rolling
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Rolling -> ShowS
showsPrec :: Int -> Rolling -> ShowS
$cshow :: Rolling -> String
show :: Rolling -> String
$cshowList :: [Rolling] -> ShowS
showList :: [Rolling] -> ShowS
Show)

-- | Full trailing windows only; variance is population variance (ddof=0).
rolling :: Rolling -> Int -> [Double] -> Either MathError [Double]
rolling :: Rolling -> Int -> [Double] -> Either MathError [Double]
rolling Rolling
method Int
width [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    Bool -> String -> Either MathError ()
require
        (Int
width Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
width Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Bool -> Bool -> Bool
&& Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
width Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
4000000)
        String
"invalid rolling window/work limit"
    [Double] -> Either MathError [Double]
checked [[Double] -> Double
apply (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
width (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
i [Double]
xs)) | Int
i <- [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
width]]
  where
    apply :: [Double] -> Double
apply [Double]
ys = case Rolling
method of
        Rolling
MovingMean -> [Double] -> Double
average [Double]
ys
        Rolling
MovingVariance -> [Double] -> Double
var [Double]
ys
        Rolling
MovingSD -> Double -> Double
forall a. Floating a => a -> a
sqrt ([Double] -> Double
var [Double]
ys)
        Rolling
MovingMinimum -> [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
ys
        Rolling
MovingMaximum -> [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
ys
    var :: [Double] -> Double
var [Double]
ys = let m :: Double
m = [Double] -> Double
average [Double]
ys in [Double] -> Double
average [(Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
m) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) | Double
x <- [Double]
ys]

-- | Recursive EMA, initialized at the first observation; 0<alpha<=1.
exponentialAverage :: Double -> [Double] -> Either MathError [Double]
exponentialAverage :: Double -> [Double] -> Either MathError [Double]
exponentialAverage Double
alpha [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
alpha Bool -> Bool -> Bool
&& Double
alpha Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
alpha Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1) String
"EMA alpha outside (0,1]"
    [Double] -> Either MathError [Double]
checked ((Double -> Double -> Double) -> Double -> [Double] -> [Double]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\Double
a Double
x -> (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
alpha) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x) (([Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) [Double]
xs) ((Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1) [Double]
xs))

-- | Nonnegative chronological weights: first weight multiplies oldest sample.
weightedAverage :: [Double] -> [Double] -> Either MathError [Double]
weightedAverage :: [Double] -> [Double] -> Either MathError [Double]
weightedAverage [Double]
ws [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
ws
    Bool -> String -> Either MathError ()
require
        ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0) [Double]
ws Bool -> Bool -> Bool
&& Double -> Bool
finite ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
ws) Bool -> Bool -> Bool
&& [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
ws Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ws Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs)
        String
"invalid moving weights"
    Bool -> String -> Either MathError ()
require (Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ws) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
4000000) String
"weighted average work limit"
    [Double] -> Either MathError [Double]
checked
        [[Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
w Double
x -> (Double
w Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
ws) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x) [Double]
ws (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
i [Double]
xs)) | Int
i <- [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ws]]

-- | Cumulative simple return relative to the first strictly positive price.
cumulativeReturns :: [Double] -> Either MathError [Double]
cumulativeReturns :: [Double] -> Either MathError [Double]
cumulativeReturns [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    Bool -> String -> Either MathError ()
require ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) [Double]
xs) String
"positive prices required"
    [Double] -> Either MathError [Double]
checked [Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ([Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) [Double]
xs Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1 | Double
x <- [Double]
xs]

-- | Consecutive log returns, using log prices to avoid ratio overflow.
logarithmicReturns :: [Double] -> Either MathError [Double]
logarithmicReturns :: [Double] -> Either MathError [Double]
logarithmicReturns [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    Bool -> String -> Either MathError ()
require ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) [Double]
xs) String
"positive prices required"
    [Double] -> Either MathError [Double]
checked ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Floating a => a -> a
log ((Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1) [Double]
xs)) ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Floating a => a -> a
log [Double]
xs))

-- | Apply (1-B^lag)^order. Zero order preserves data; at least one output is required.
difference :: Int -> Int -> [Double] -> Either MathError [Double]
difference :: Int -> Int -> [Double] -> Either MathError [Double]
difference Int
order Int
lag [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
xs
    Bool -> String -> Either MathError ()
require
        (Int
order Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
order Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
8 Bool -> Bool -> Bool
&& Int
lag Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
order Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
lag Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs))
        String
"invalid differencing order/lag"
    [Double] -> Either MathError [Double]
checked (([Double] -> [Double]) -> [Double] -> [[Double]]
forall a. (a -> a) -> a -> [a]
iterate (\[Double]
ys -> (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
lag [Double]
ys) [Double]
ys) [Double]
xs [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
order)

-- | Biased covariance divided by lag-zero covariance, global sample mean.
autocorrelation :: Int -> [Double] -> Either MathError [Double]
autocorrelation :: Int -> [Double] -> Either MathError [Double]
autocorrelation Int
lags [Double]
xs = do
    Int -> [Double] -> Either MathError ()
validate Int
2 [Double]
xs
    Bool -> String -> Either MathError ()
require
        (Int
lags Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
lags Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Bool -> Bool -> Bool
&& Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
lags Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
4000000)
        String
"invalid ACF lags/work limit"
    let m :: Double
m = [Double] -> Double
average [Double]
xs; zs :: [Double]
zs = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Num a => a -> a -> a
subtract Double
m) [Double]
xs; den :: Double
den = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) [Double]
zs)
    Bool -> String -> Either MathError ()
require (Double
den Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
den) String
"ACF requires positive finite variance"
    [Double] -> Either MathError [Double]
checked [[Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
k [Double]
zs) [Double]
zs) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
den | Int
k <- [Int
0 .. Int
lags]]

-- | Corr(x_t,y_(t+lag)), global means and full lag-zero energy normalization.
crossCorrelation :: Int -> [Double] -> [Double] -> Either MathError Double
crossCorrelation :: Int -> [Double] -> [Double] -> Either MathError Double
crossCorrelation Int
lag [Double]
xs [Double]
ys = do
    Int -> [Double] -> Either MathError ()
validate Int
2 [Double]
xs
    Int -> [Double] -> Either MathError ()
validate Int
2 [Double]
ys
    Bool -> String -> Either MathError ()
require
        ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys Bool -> Bool -> Bool
&& Integer -> Integer
forall a. Num a => a -> a
abs (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
lag) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs))
        String
"cross-correlation shape/lag"
    let a :: [Double]
a = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Num a => a -> a -> a
subtract ([Double] -> Double
average [Double]
xs)) [Double]
xs
        b :: [Double]
b = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Num a => a -> a -> a
subtract ([Double] -> Double
average [Double]
ys)) [Double]
ys
        energy :: [a] -> a
energy [a]
values = a -> a
forall a. Floating a => a -> a
sqrt ([a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) [a]
values))
        den :: Double
den = [Double] -> Double
forall {a}. Floating a => [a] -> a
energy [Double]
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall {a}. Floating a => [a] -> a
energy [Double]
b
        ([Double]
u, [Double]
v) = if Int
lag Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 then ([Double]
a, Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
lag [Double]
b) else (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop (-Int
lag) [Double]
a, [Double]
b)
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
den Bool -> Bool -> Bool
&& Double
den Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"cross-correlation requires finite positive variances"
    ([Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) ([Double] -> Double)
-> Either MathError [Double] -> Either MathError Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [Double] -> Either MathError [Double]
checked [[Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
u [Double]
v) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
den]

{- | Levinson recursion for a positive definite Toeplitz covariance. Returns
AR coefficients, innovation variance, and PACF (without the lag-zero one).
-}
durbinLevinson :: [Double] -> Either MathError ([Double], Double, [Double])
durbinLevinson :: [Double] -> Either MathError ([Double], Double, [Double])
durbinLevinson [Double]
covariance = do
    Int -> [Double] -> Either MathError ()
validate Int
1 [Double]
covariance
    Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
covariance Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
257 Bool -> Bool -> Bool
&& ([Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) [Double]
covariance Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"invalid covariance/order"
    Int
-> [Double]
-> Double
-> [Double]
-> Either MathError ([Double], Double, [Double])
go Int
1 [] (([Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) [Double]
covariance) []
  where
    go :: Int
-> [Double]
-> Double
-> [Double]
-> Either MathError ([Double], Double, [Double])
go Int
k [Double]
phi Double
variance [Double]
pacf
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
covariance = ([Double], Double, [Double])
-> Either MathError ([Double], Double, [Double])
forall a b. b -> Either a b
Right ([Double]
phi, Double
variance, [Double]
pacf)
        | Bool
otherwise = do
            let reflection :: Double
reflection = ([Double]
covariance [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[Double]
phi [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
covariance [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
j) | Int
j <- [Int
1 .. Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
variance
                next :: [Double]
next = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
a Double
b -> Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
reflection Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b) [Double]
phi ([Double] -> [Double]
forall a. [a] -> [a]
reverse [Double]
phi) [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double
reflection]
                v :: Double
v = Double
variance Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
reflection Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
reflection)
            _ <- [Double] -> Either MathError [Double]
checked (Double
reflection Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: Double
v Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [Double]
next)
            require (abs reflection < 1 && v > 0) "Toeplitz covariance is not positive definite"
            go (k + 1) next v (pacf ++ [reflection])

-- | Biased sample autocovariances with centered data; returns mean separately.
yuleWalker :: Int -> [Double] -> Either MathError ([Double], Double, Double)
yuleWalker :: Int -> [Double] -> Either MathError ([Double], Double, Double)
yuleWalker Int
order [Double]
xs = do
    acf <- Int -> [Double] -> Either MathError [Double]
autocorrelation Int
order [Double]
xs
    let m = [Double] -> Double
average [Double]
xs; variance = [Double] -> Double
average [(Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
m) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) | Double
x <- [Double]
xs]
    (phi, v, _) <- durbinLevinson (map (* variance) acf)
    Right (phi, v, m)

{- | Given maximum lag then samples, compute ACF and Durbin-Levinson partial correlations,
returning lag zero (1) followed by requested PACF lags. Invalid lag/data or singular
covariance recurrence returns Left MathError.
-}
partialAutocorrelation :: Int -> [Double] -> Either MathError [Double]
partialAutocorrelation :: Int -> [Double] -> Either MathError [Double]
partialAutocorrelation Int
order [Double]
xs = do
    acf <- Int -> [Double] -> Either MathError [Double]
autocorrelation Int
order [Double]
xs
    (_, _, pacf) <- durbinLevinson acf
    Right (1 : pacf)

-- Relative position between finite ordered anchors, avoiding opposite-sign
-- endpoint subtraction overflow.
fraction :: Double -> Double -> Double -> Double
fraction :: Double -> Double -> Double -> Double
fraction Double
a Double
b Double
t = if Double -> Bool
finite (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) then (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) else (Double
t Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
b Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)