plans.stochastic.arma#

Trend, cycle and AR(1) memory model for synthetic series.

A series is the sum of a linear trend, a mean-zero cycle and an AR(1) remainder driven by Gaussian shocks. For steps \(t = 0, 1, \dots, n-1\):

\[x_t = \underbrace{b + a t}_{\text{trend}} + \underbrace{A (g_t - m_k)}_{\text{cycle}} + R_t, \qquad R_t = \varphi R_{t-1} + \varepsilon_t, \qquad \varepsilon_t \sim \mathcal{N}(0, \sigma^2)\]

Parameters map to symbols as start \(b\), slope \(a\), amplitude \(A\), period \(P\), phase \(t_0\), power \(k\), phi \(\varphi\) and sigma \(\sigma\). Functions take plain arrays and a period (steps per cycle), so they know nothing about dates or units.

Three operations link observed data, parameters and synthetic data: fit_params() goes from data to parameters, decompose() splits data into components given the parameters, and build() makes a series from parameters.

Functions

`ar1_from_shocks`(shocks, phi)	Get the AR(1) remainder from given shocks, with the first value set to the first shock.
`average_from_start`(start, slope, n)	Get the trend average over the record from the trend start.
`build`(params, n[, seed, shocks])	Make a series from the parameters.
`cycle_component`(n, amplitude, period, phase, ...)	Get the cycle term, which averages zero over a cycle.
`cycle_shape`(n, period, phase, power)	Get the shape of the cycle, a raised cosine to a power.
`decompose`(x, params)	Split a series into components, given the parameters.
`describe`(params, n[, positive])	Get a printout with all views of the same model.
`draw_shocks`(n, sigma, rng)	Draw independent Gaussian shocks.
`fit_params`(x, period)	Fit the model to a series.
`get_floor`(params, n)	Get the lowest value of trend plus cycle, \(\min_t (T_t + C_t)\).
`make_params`(period[, start, average, slope, ...])	Make the parameters from what is known, converting to the stored form.
`shape_mean`(power)	Get the average of the shape over a cycle.
`simulate_ar1`(n, phi, sigma, rng)	Simulate an AR(1) remainder that is stationary from the first step.
`slope_from_cycle_change`(change, period)	Get the slope from the trend change per cycle \(\Delta T\).
`slope_from_end`(start, end, n)	Get the slope from the trend values at the first and last steps.
`start_from_average`(average, slope, n)	Get the trend start from the trend average over the record.
`trend_line`(n, start, slope)	Get the linear trend.

plans.stochastic.arma.trend_line(n, start, slope)[source]#

Get the linear trend.

\[T_t = b + a t\]

Parameters:

n (int) – number of steps
start (float) – trend value at step 0
slope (float) – trend change per step

Returns:

trend values

Return type:

numpy.ndarray

plans.stochastic.arma.cycle_shape(n, period, phase, power)[source]#

Get the shape of the cycle, a raised cosine to a power.

\[g_t = \left[ \frac{1 + \cos\theta_t}{2} \right]^{k}, \qquad \theta_t = \frac{2 \pi (t - t_0)}{P}\]

Parameters:

n (int) – number of steps
period (float) – steps per cycle
phase (float) – step at which the cycle peaks
power (float) – shape. 1 gives a plain cosine, larger gives a sharper peak.

Returns:

shape values, from 0 to 1

Return type:

numpy.ndarray

plans.stochastic.arma.shape_mean(power)[source]#

Get the average of the shape over a cycle.

\[m_k = \frac{\Gamma(k + \tfrac{1}{2})}{\sqrt{\pi} \, \Gamma(k + 1)}\]

It equals 0.5 for \(k = 1\) and 0.375 for \(k = 2\).

Parameters:: power (float) – shape
Returns:: average of the shape
Return type:: float

plans.stochastic.arma.cycle_component(n, amplitude, period, phase, power)[source]#

Get the cycle term, which averages zero over a cycle.

\[C_t = A (g_t - m_k)\]

For \(k = 1\) this is a plain cosine of half-range \(A / 2\).

Parameters:

n (int) – number of steps
amplitude (float) – peak-to-trough range of the cycle term
period (float) – steps per cycle
phase (float) – step at which the cycle peaks
power (float) – shape

Returns:

cycle values

Return type:

numpy.ndarray

plans.stochastic.arma.draw_shocks(n, sigma, rng)[source]#

Draw independent Gaussian shocks.

\[\varepsilon_t \sim \mathcal{N}(0, \sigma^2)\]

Parameters:

n (int) – number of steps
sigma (float) – shock standard deviation
rng (numpy.random.Generator) – random generator

Returns:

shocks

Return type:

numpy.ndarray

plans.stochastic.arma.ar1_from_shocks(shocks, phi)[source]#

Get the AR(1) remainder from given shocks, with the first value set to the first shock.

\[R_0 = \varepsilon_0, \qquad R_t = \varphi R_{t-1} + \varepsilon_t\]

Parameters:

shocks (numpy.ndarray) – shocks
phi (float) – AR(1) memory

Returns:

remainder

Return type:

numpy.ndarray

plans.stochastic.arma.simulate_ar1(n, phi, sigma, rng)[source]#

Simulate an AR(1) remainder that is stationary from the first step.

The first shock is scaled to the stationary variance instead of discarding a burn-in. Draws therefore stay aligned when phi changes, so the same generator state isolates the effect of phi.

\[R_0 = \frac{\varepsilon_0}{\sqrt{1 - \varphi^2}}, \qquad \mathrm{Var}(R_t) = \frac{\sigma^2}{1 - \varphi^2}\]

Parameters:

n (int) – number of steps
phi (float) – AR(1) memory, inside (-1, 1)
sigma (float) – shock standard deviation
rng (numpy.random.Generator) – random generator

Returns:

remainder

Return type:

numpy.ndarray

plans.stochastic.arma.fit_params(x, period)[source]#

Fit the model to a series.

Trend and cycle are fitted jointly by nonlinear least squares from several starting shapes, keeping the best. power is bounded by POWER_MIN and POWER_MAX, so a cycle with a sharp trough and a flat top cannot be represented.

Memory and shocks come from the remainder \(R_t\), after removing its mean. The memory is a regression through the origin, and \(\hat{\sigma}\) is the standard deviation of the residuals of that regression:

\[\hat{\varphi} = \frac{\sum_t R_{t-1} R_t}{\sum_t R_{t-1}^2}, \qquad \hat{\sigma} = \mathrm{sd}\left( R_t - \hat{\varphi} R_{t-1} \right)\]

Warning

Trend and cycle estimates ignore the memory, so they are less precise than they could be, and phi tends to be low on short series.

Parameters:

x (numpy.ndarray) – series without nan values
period (float) – steps per cycle

Returns:

parameters. Keys: start, slope, amplitude, period, phase, power, phi, sigma

Return type:

dict

plans.stochastic.arma.decompose(x, params)[source]#

Split a series into components, given the parameters.

\[R_t = x_t - T_t - C_t, \qquad \varepsilon_0 = R_0, \qquad \varepsilon_t = R_t - \varphi R_{t-1}\]

Parameters:

x (numpy.ndarray) – series
params (dict) – parameters. See fit_params().

Returns:

arrays by component. Keys: trend, cycle, remainder, shocks

Return type:

dict

plans.stochastic.arma.build(params, n, seed=None, shocks=None)[source]#

Make a series from the parameters.

To isolate the effect of one parameter, change it and build again with the same seed or the same shocks.

Parameters:

params (dict) – parameters. See fit_params().
n (int) – number of steps
seed (int or None) – random seed, used when shocks is None
shocks (numpy.ndarray or None) – shocks to use instead of drawing them, of size n

Returns:

series

Return type:

numpy.ndarray

plans.stochastic.arma.get_floor(params, n)[source]#

Get the lowest value of trend plus cycle, \(\min_t (T_t + C_t)\).

Note

This is the deterministic floor. Shocks are Gaussian and additive, so a synthetic series can still fall below it.

Parameters:

params (dict) – parameters. See fit_params().
n (int) – number of steps

Returns:

lowest value of trend plus cycle

Return type:

float

plans.stochastic.arma.describe(params, n, positive=False)[source]#

Get a printout with all views of the same model.

Baseline and topline are the extremes of the deterministic cycle, at the trend average \(\bar{T}\). Shocks move the observed extremes beyond them.

\[\text{baseline} = \bar{T} - A m_k, \qquad \text{topline} = \bar{T} + A (1 - m_k)\]

Parameters:

params (dict) – parameters. See fit_params().
n (int) – number of steps, used by the views that depend on the length
positive (bool) – set True for variables that must stay positive. It warns if the deterministic floor is below zero.

Returns:

formatted text

Return type:

str

plans.stochastic.arma.make_params(period, start=None, average=None, slope=None, end=None, cycle_change=None, n=None, amplitude=None, baseline=None, topline=None, phase=0.0, power=1.0, phi=0.0, sigma=0.0)[source]#

Make the parameters from what is known, converting to the stored form.

The trend level comes from start or average, its change from slope, end or cycle_change, and the cycle from amplitude or from both baseline and topline, with \(A = \text{topline} - \text{baseline}\). Giving more than one option for the same item raises an error. Length-based conversions apply once here, so a longer series continues the same line.

See start_from_average(), slope_from_end() and slope_from_cycle_change() for the formulas.

Parameters:

period (float) – steps per cycle
start (float) – trend value at step 0
average (float) – trend average over the record. Requires n.
slope (float) – trend change per step
end (float) – trend value at the last step. Requires start and n.
cycle_change (float) – trend change per cycle
n (int) – number of steps, used only by the length-based conversions
amplitude (float) – peak-to-trough range of the cycle term
baseline (float) – lowest value of the deterministic cycle
topline (float) – highest value of the deterministic cycle
phase (float) – step at which the cycle peaks
power (float) – shape
phi (float) – AR(1) memory, inside (-1, 1)
sigma (float) – shock standard deviation

Returns:

parameters. See fit_params().

Return type:

dict

plans.stochastic.arma.start_from_average(average, slope, n)[source]#

Get the trend start from the trend average over the record.

\[b = \bar{T} - a \frac{n - 1}{2}\]

Parameters:

average (float) – trend average
slope (float) – trend change per step
n (int) – number of steps

Returns:

trend value at step 0

Return type:

float

plans.stochastic.arma.average_from_start(start, slope, n)[source]#

Get the trend average over the record from the trend start.

\[\bar{T} = b + a \frac{n - 1}{2}\]

Parameters:

start (float) – trend value at step 0
slope (float) – trend change per step
n (int) – number of steps

Returns:

trend average

Return type:

float

plans.stochastic.arma.slope_from_end(start, end, n)[source]#

Get the slope from the trend values at the first and last steps.

\[a = \frac{T_{n-1} - b}{n - 1}\]

Parameters:

start (float) – trend value at step 0
end (float) – trend value at the last step
n (int) – number of steps

Returns:

trend change per step

Return type:

float

plans.stochastic.arma.slope_from_cycle_change(change, period)[source]#

Get the slope from the trend change per cycle \(\Delta T\).

\[a = \frac{\Delta T}{P}\]

Parameters:

change (float) – trend change per cycle
period (float) – steps per cycle

Returns:

trend change per step

Return type:

float