Levinson-Durbin and Autocovariance

As discussed, I am implementing estimation methods for `phi` in `AR` modelling. We covered `yule_walker` earlier, I’ll write a post about that. After it’s implementation, we go ahead with another estimation method - Levinson Durbin

Levinson-Durbin requires timeline series to be demeaned(`series = series - series.mean`) and it’s autocovavirance.

Autocovariance of series is represented by summation of summation of product of series with series at lag `k`. That is, summation of `(x_i * x_{i+lag})`. It is also directly related with `acf` of series as `acf(k) = acvf(h) / acvf(0)`. It’s code can now be found in `Statsample::TimeSeries`’s `acvf` method.

Now, with the help of autocovariance series, our `levinson_durbin` function recursively computes the following parameters:

• sigma_v : estimation of error variance
• arcoefs : AR phi values for timeseries
• pac : unbiased levinson pacf estiation
• sigma : sigma for AR.

L-D performs recursive matrix and vector multiplications to populate it’s toeplitz matrix. Here is some code depicting those manipulations:

Implementation can be found here.

Now, in this week, I will integrate this in AR modelling and perform some tests to verify the estimation. And will soon start with next estimation method :)

Cheers,
Ankur Goel

AR/MA, ARMA Acf - Pacf Visualizations

As mentioned in previous post, I have been working with Autoregressive and Moving Average simulations.
To test the correctness of estimations by our simulations, we employ acf(Autocorrelation) and pacf(partial autocorrelation) to our use. For different order of AR and MA, we get the varying visualizations with them, such as:

• Exponential decreasing curves.
• Damped sine waves.
• Positive and negative spikes, etc.

While analyzing and writing tests for same, I also took some time to visualize that data on ilne and bar charts to get a clearer picture:

AR(1) process

AR(1) process is the autoregressive simulation with order p = 1, i.e, with one value of phi.
Ideal AR(p) process is represented by:

To simulate this, install `statsample-timeseries` from here.

Here, number of observations, n = 1500 (greater value is preferrable for best fit), p = 1, with phi = [0.9].

ACF

To generate it’s autocorrelation

For an AR(1) process, `acf` must exponentially decay if `phi > 0`, or alternate in sign if `phi < 0`Ref. Go through the analysis above. It can be visualized as:
When `phi > 0`, acf decreases exponentially:

When `phi < 0`, you get the alternate acf lags:

PACF

To generate it’s partial autocorrelation:

For AR(1) process, `pacf` must have a spike at lag 1, then 0. Former spike must be positive if `phi > 0`, otherwise, negative spike. Have a look at the pacf series generated above. On visualizing the data:
When `phi > 0`, positive lag at 1 and 0(contains 1.0):

When `phi < 0`, negative lag at 1:

Here is the representation of ideal acf-vs-pacf for positive phi in AR(1):

AR(P) process

Simulation of AR(p) process is similar as AR(1).

ACF

For AR(p), `acf` must give a damping sine wave. The pattern is greatly dependent on the value and sign of phi parameters.
When positive content in phi coefficients is more, you will get a sine wave starting from positive side, else, sine wave will start from negative side.
Notice, the damping sine wave starting from positive side here:

and negative side here..

PACF

`pacf` gives spike at lag 0(value = 1.0, default) and from lag 1 to lag k. The example above, features AR(2) process, for this, we must get spikes at lag 1 - 2 as:

MA(1) process

MA(1) process is the moving average simulation with order `q = 1`, i.e, with one value of `theta`.
To simulate this, use `ma_sim` method from `Statsample::ARIMA::ARIMA`

For `theta > 0`, for `MA(1)`, we must get a positive spike at lag 1 as:
For `theta < 0`, the spike at lag 1 must be in negatie direction as:

When I put these two visualizations aside each other, the visualization seems quite fit:

MA(q) process

MA(q) process. `Order = q` => Number of theta coefficients = q.
Ideal MA(q) process is represented by:

ACF

Similar to AR(1) simulation, it will have spikes for lag 1 - lag p as :

PACF

In pacf of MA(q) simulation, we observe exponentially decaying/damping sine wave.

ARMA(p, q) process

ARMA(p, q) is combination of autoregressive and moving average simulations.
When `q = 0`, the process is called as pure autoregressive process; when `p = 0`, the process is purely moving average.
The simulator of ARMA can be found as `arma_sim` in `Statsample::ARIMA::ARIMA`.
For `ARMA(1, 1)` process, here are the comparisons of the visualizations from `R` and this code, which just made my day :)

Quite Fit!

Cheers,
- Ankur Goel

AR and MA Simulations

I have been reading and coding AR(Autoregressive Model) and Moving Average(MA) basics.
First, I am very grateful for all the patience by Claudio Bustos to help me with understanding this. There is still lot to be learnt in these models, but I will keep bugging him. :D

We wrote the simulations for `AR(p)` and `MA(q)`. The idea behind creating them is to first simulate the process with known coefficients and then move on to write ARMA process which could also estimate the coefficients.

Autoregressive Model

This model is represented by AR(p) where p is the order of this model. For a pure autoregressive model, we consider order of moving average model as 0.

`AR(p)` is represented by the following equation:

Here, are the parameters of model, and is the error noise.
To realize this model, we have to observe and keep track of previous `x(t)` values, and realize current value with observed summation and error noise.

Here is that code fragment of general AR(p) simulator:

• It creates a buffer initially to prevent ‘broken values’, and trims that buffer before returning the result.
• It goes on to create the backshifts vector. The backshift vector depends on the ongoing iteration from `(1..n)`.
• After creating backshifts vector, it extracts the required `phi` values from the stack.
• Then it performs vector multiplication - `backshifts * parameters`, then adds the result.
• The current `x(t)` values is then returned with active white noise.

There are certain tests which will now be performed in context of `acf` and `pacf`, which I earlier coded. They form the basis of correctness of our autoregressive estimation. :) Expect the next post about those tests.

Moving Average Model

This model is represented by MA(q), where `q` is the order of this model. Again, for pure moving-average model, we consider `p=0`.

This is represented by following equation:

Unlike autoregressive model, this model was somewhat hard to obtain. It needs to obsere previous error noise instead of previous `x(t)` values. And the series largely depends upon the order of the model.

It’s code can also be found on my github branch.

I am now currently working on those tests analysis, I mentioned; for various combinations of `AR(p)` and `MA(q)`. As they get over, I will move to realize few algorithms like yule walker and burgs for estimation of coefficients.

Cheers,
-Ankur Goel