Periodic Levinson-Durbin algorithm
pcalg1.RdCalculate partial periodic autocorrelations, forward and backward prediction coefficients and error variances using the periodic Levinson-Durbin algorithm.
Arguments
- r
periodic autocovariances, a matrix, see `Details'.
- p
autoregressive orders, numeric vector.
Details
alg1(r,p) calculates the partial periodic correlations from
autocovariances r and autoregression orders p. The
matrix r has the same format as that of the r slot of
pcAcvf objects. The periodicity, d, is set equal to
the number of rows in r. If the length of p is not
equal to the periodicity, all autoregressive orders are set to the
first element of p. This last feature is really meant to be
used only with a scalar p.
The convention for the signs of the coefficients is the one from Boshnakov(1996) and is consistent with other R time series functions.
pmax below stands for the maximal element of p,
i.e. the maximal AR order.
As in the non-periodic case, the periodic Levinson-Durbin algorithm
fits recursively models of order 0, 1, ..., pmax. Namely, at
step i the AR orders for all seasons are set to
i. This is done in a way that correctly handles the case when
not all elements of p are equal, see the references.
The essential quantities calculated by the periodic Levinson-Durbin algorithm are returned as matrices, whose \(i\)th rows contain values for season \(i\). The complete details depend on the quantities, as described below.
The partial autocorrelations, the forward innovation variances and
the backward innovation variances are returned as matrices with
d rows and 1+pmax columns, whose j-th columns contain
the quantities for order j-1 (partial autocorrelations, forward
innovation variances and backward innovation variances, respectively).
Note that the lag-0 partial autocorrelations are the autocovariances
for lag 0, see the references for details.
The forward autoregression parameters are returned as a list whose \(j\)th element is a matrix containing the coefficients for order \(j\). Similarly for the backward autoregression parameters.
One often is interested in the model of order p only. Its
coefficients are given by af[[pmax]], while the innovation
variances are in the last column of fv.
Value
A list with the following elements.
- orders
autoregression orders
- be
partial autocorrelations, a matrix with d rows
- fv
forward innovation variances, a matrix with d rows
- bv
backward innovation variances, a matrix with d rows
- af
forward autoregression parameters, a list with one element for the parameters for each order.
- ab
backward autoregression parameters, a list with one element for the parameters for each order.
References
Boshnakov GN (1996). “Recursive computation of the parameters of periodic autoregressive moving-average processes.” J. Time Ser. Anal., 17(4), 333--349. ISSN 0143-9782, doi: 10.1111/j.1467-9892.1996.tb00281.x .
Lambert-Lacroix S (2000). “On periodic autoregressive process estimation .” IEEE Transactions on Signal Processing, 48( 6 ), 1800-1803.
Lambert-Lacroix S (2005). “ Extension of autocovariance coefficients sequence for periodically correlated processes.” Journal of Time Series Analysis, 26(3), 423-435.
Note
The autoregression orders of the output are not necessarilly the same as those specified in the call. There may be no PAR model with the requested orders, see the references.
Examples
r1 <- rbind(c(1,0.81,0.729),c(1,0.90,0.900))
alg1(r1,2)
#> $orders
#> [1] 2 2
#>
#> $be
#> [,1] [,2] [,3]
#> [1,] 1 0.81 -4.343275e-16
#> [2,] 1 0.90 6.689647e-01
#>
#> $fv
#> [,1] [,2] [,3]
#> [1,] 1 0.3439 0.3439000
#> [2,] 1 0.1900 0.1049724
#>
#> $bv
#> [,1] [,2] [,3]
#> [1,] 1 0.1900 0.19
#> [2,] 1 0.3439 0.19
#>
#> $af
#> $af[[1]]
#> [,1]
#> [1,] -0.81
#> [2,] -0.90
#>
#> $af[[2]]
#> [,1] [,2]
#> [1,] -0.8100000 5.843279e-16
#> [2,] -0.4972376 -4.972376e-01
#>
#>
#> $ab
#> $ab[[1]]
#> [,1]
#> [1,] -0.90
#> [2,] -0.81
#>
#> $ab[[2]]
#> [,1] [,2]
#> [1,] -0.9 3.228331e-16
#> [2,] 0.0 -9.000000e-01
#>
#>
## pc2 <- pcAcvf(init=r1)
## pc2a <- pcAcvf(init=r1,seasonnames=c("am","pm"), periodunit="day")
# example of Lambert-Lacroix
data(ex1f)
pc3 <- slMatrix(period=2,maxlag=5,f=ex1f,type="tt")
res0p2 <- alg1(pc3[],c(0,2))
res1p2 <- alg1(pc3[],c(1,2))
res3p3 <- alg1(pc3[],c(3,3))
paramsys1 <- pcarma_param_system(pc3, NULL, NULL, 2, 0, 2)
t1 <- solve(paramsys1$A,paramsys1$b)
# this is from tests.r but I have lost t1
# set it to pc3 below
# note: t1 is not the t1 computed above and in other examples!
t1 <- pc3
t1
#> An object of class "slMatrix"
#> Slot "m":
#> lag
#> season 0 1 2 3 4 5
#> 1 1 0.81 0.729 0.729 0.6561 0.6561
#> 2 1 0.90 0.900 0.810 0.8100 0.7290
#>
t1[]
#> lag
#> season 0 1 2 3 4 5
#> 1 1 0.81 0.729 0.729 0.6561 0.6561
#> 2 1 0.90 0.900 0.810 0.8100 0.7290
alg1(t1[],c(1,1))
#> $orders
#> [1] 1 1
#>
#> $be
#> [,1] [,2]
#> [1,] 1 0.81
#> [2,] 1 0.90
#>
#> $fv
#> [,1] [,2]
#> [1,] 1 0.3439
#> [2,] 1 0.1900
#>
#> $bv
#> [,1] [,2]
#> [1,] 1 0.1900
#> [2,] 1 0.3439
#>
#> $af
#> $af[[1]]
#> [,1]
#> [1,] -0.81
#> [2,] -0.90
#>
#>
#> $ab
#> $ab[[1]]
#> [,1]
#> [1,] -0.90
#> [2,] -0.81
#>
#>
alg1(t1[],c(1,0))
#> $orders
#> [1] 1 0
#>
#> $be
#> [,1] [,2]
#> [1,] 1 0.81
#> [2,] 1 0.00
#>
#> $fv
#> [,1] [,2]
#> [1,] 1 0.3439
#> [2,] 1 1.0000
#>
#> $bv
#> [,1] [,2]
#> [1,] 1 1.0000
#> [2,] 1 0.3439
#>
#> $af
#> $af[[1]]
#> [,1]
#> [1,] -0.81
#> [2,] 0.00
#>
#>
#> $ab
#> $ab[[1]]
#> [,1]
#> [1,] 0.00
#> [2,] -0.81
#>
#>
alg1(t1[],c(0,1))
#> $orders
#> [1] 0 1
#>
#> $be
#> [,1] [,2]
#> [1,] 1 0.0
#> [2,] 1 0.9
#>
#> $fv
#> [,1] [,2]
#> [1,] 1 1.00
#> [2,] 1 0.19
#>
#> $bv
#> [,1] [,2]
#> [1,] 1 0.19
#> [2,] 1 1.00
#>
#> $af
#> $af[[1]]
#> [,1]
#> [1,] 0.0
#> [2,] -0.9
#>
#>
#> $ab
#> $ab[[1]]
#> [,1]
#> [1,] -0.9
#> [2,] 0.0
#>
#>
alg1(t1[],c(5,5))
#> $orders
#> [1] 5 5
#>
#> $be
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 1 0.81 0.0000000 0 -4.343275e-16 -4.343275e-16
#> [2,] 1 0.90 0.6689647 0 0.000000e+00 7.861328e-16
#>
#> $fv
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 1 0.3439 0.3439000 0.3439000 0.3439000 0.3439000
#> [2,] 1 0.1900 0.1049724 0.1049724 0.1049724 0.1049724
#>
#> $bv
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 1 0.1900 0.19 0.19 0.19 0.19
#> [2,] 1 0.3439 0.19 0.19 0.19 0.19
#>
#> $af
#> $af[[1]]
#> [,1]
#> [1,] -0.81
#> [2,] -0.90
#>
#> $af[[2]]
#> [,1] [,2]
#> [1,] -0.8100000 0.0000000
#> [2,] -0.4972376 -0.4972376
#>
#> $af[[3]]
#> [,1] [,2] [,3]
#> [1,] -0.8100000 0.0000000 0
#> [2,] -0.4972376 -0.4972376 0
#>
#> $af[[4]]
#> [,1] [,2] [,3] [,4]
#> [1,] -0.8100000 0.0000000 -5.258951e-16 5.843279e-16
#> [2,] -0.4972376 -0.4972376 0.000000e+00 0.000000e+00
#>
#> $af[[5]]
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] -0.8100000 0.0000000 -1.05179e-15 5.843279e-16 5.843279e-16
#> [2,] -0.4972376 -0.4972376 0.00000e+00 5.258951e-16 -5.843279e-16
#>
#>
#> $ab
#> $ab[[1]]
#> [,1]
#> [1,] -0.90
#> [2,] -0.81
#>
#> $ab[[2]]
#> [,1] [,2]
#> [1,] -0.9 0.0
#> [2,] 0.0 -0.9
#>
#> $ab[[3]]
#> [,1] [,2] [,3]
#> [1,] -0.9 0.0 0
#> [2,] 0.0 -0.9 0
#>
#> $ab[[4]]
#> [,1] [,2] [,3] [,4]
#> [1,] -0.9 0.0 -2.614948e-16 3.228331e-16
#> [2,] 0.0 -0.9 0.000000e+00 0.000000e+00
#>
#> $ab[[5]]
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] -9.000000e-01 0.0 2.644003e-16 8.487282e-16 -1.057634e-15
#> [2,] 1.886404e-31 -0.9 0.000000e+00 -2.614948e-16 3.228331e-16
#>
#>
alg1(t1[],c(2,2))
#> $orders
#> [1] 2 2
#>
#> $be
#> [,1] [,2] [,3]
#> [1,] 1 0.81 0.0000000
#> [2,] 1 0.90 0.6689647
#>
#> $fv
#> [,1] [,2] [,3]
#> [1,] 1 0.3439 0.3439000
#> [2,] 1 0.1900 0.1049724
#>
#> $bv
#> [,1] [,2] [,3]
#> [1,] 1 0.1900 0.19
#> [2,] 1 0.3439 0.19
#>
#> $af
#> $af[[1]]
#> [,1]
#> [1,] -0.81
#> [2,] -0.90
#>
#> $af[[2]]
#> [,1] [,2]
#> [1,] -0.8100000 0.0000000
#> [2,] -0.4972376 -0.4972376
#>
#>
#> $ab
#> $ab[[1]]
#> [,1]
#> [1,] -0.90
#> [2,] -0.81
#>
#> $ab[[2]]
#> [,1] [,2]
#> [1,] -0.9 0.0
#> [2,] 0.0 -0.9
#>
#>
alg1(t1[],c(2,3))
#> $orders
#> [1] 2 3
#>
#> $be
#> [,1] [,2] [,3] [,4]
#> [1,] 1 0.81 0.0000000 0
#> [2,] 1 0.90 0.6689647 0
#>
#> $fv
#> [,1] [,2] [,3] [,4]
#> [1,] 1 0.3439 0.3439000 0.3439000
#> [2,] 1 0.1900 0.1049724 0.1049724
#>
#> $bv
#> [,1] [,2] [,3] [,4]
#> [1,] 1 0.1900 0.19 0.19
#> [2,] 1 0.3439 0.19 0.19
#>
#> $af
#> $af[[1]]
#> [,1]
#> [1,] -0.81
#> [2,] -0.90
#>
#> $af[[2]]
#> [,1] [,2]
#> [1,] -0.8100000 0.0000000
#> [2,] -0.4972376 -0.4972376
#>
#> $af[[3]]
#> [,1] [,2] [,3]
#> [1,] -0.8100000 0.0000000 0
#> [2,] -0.4972376 -0.4972376 0
#>
#>
#> $ab
#> $ab[[1]]
#> [,1]
#> [1,] -0.90
#> [2,] -0.81
#>
#> $ab[[2]]
#> [,1] [,2]
#> [1,] -0.9 0.0
#> [2,] 0.0 -0.9
#>
#> $ab[[3]]
#> [,1] [,2] [,3]
#> [1,] -0.9 0.0 0
#> [2,] 0.0 -0.9 0
#>
#>
alg1(t1[],c(3,3))
#> $orders
#> [1] 3 3
#>
#> $be
#> [,1] [,2] [,3] [,4]
#> [1,] 1 0.81 0.0000000 0
#> [2,] 1 0.90 0.6689647 0
#>
#> $fv
#> [,1] [,2] [,3] [,4]
#> [1,] 1 0.3439 0.3439000 0.3439000
#> [2,] 1 0.1900 0.1049724 0.1049724
#>
#> $bv
#> [,1] [,2] [,3] [,4]
#> [1,] 1 0.1900 0.19 0.19
#> [2,] 1 0.3439 0.19 0.19
#>
#> $af
#> $af[[1]]
#> [,1]
#> [1,] -0.81
#> [2,] -0.90
#>
#> $af[[2]]
#> [,1] [,2]
#> [1,] -0.8100000 0.0000000
#> [2,] -0.4972376 -0.4972376
#>
#> $af[[3]]
#> [,1] [,2] [,3]
#> [1,] -0.8100000 0.0000000 0
#> [2,] -0.4972376 -0.4972376 0
#>
#>
#> $ab
#> $ab[[1]]
#> [,1]
#> [1,] -0.90
#> [2,] -0.81
#>
#> $ab[[2]]
#> [,1] [,2]
#> [1,] -0.9 0.0
#> [2,] 0.0 -0.9
#>
#> $ab[[3]]
#> [,1] [,2] [,3]
#> [1,] -0.9 0.0 0
#> [2,] 0.0 -0.9 0
#>
#>
alg1(t1[],c(4,4))
#> $orders
#> [1] 4 4
#>
#> $be
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1 0.81 0.0000000 0 -4.343275e-16
#> [2,] 1 0.90 0.6689647 0 0.000000e+00
#>
#> $fv
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1 0.3439 0.3439000 0.3439000 0.3439000
#> [2,] 1 0.1900 0.1049724 0.1049724 0.1049724
#>
#> $bv
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1 0.1900 0.19 0.19 0.19
#> [2,] 1 0.3439 0.19 0.19 0.19
#>
#> $af
#> $af[[1]]
#> [,1]
#> [1,] -0.81
#> [2,] -0.90
#>
#> $af[[2]]
#> [,1] [,2]
#> [1,] -0.8100000 0.0000000
#> [2,] -0.4972376 -0.4972376
#>
#> $af[[3]]
#> [,1] [,2] [,3]
#> [1,] -0.8100000 0.0000000 0
#> [2,] -0.4972376 -0.4972376 0
#>
#> $af[[4]]
#> [,1] [,2] [,3] [,4]
#> [1,] -0.8100000 0.0000000 -5.258951e-16 5.843279e-16
#> [2,] -0.4972376 -0.4972376 0.000000e+00 0.000000e+00
#>
#>
#> $ab
#> $ab[[1]]
#> [,1]
#> [1,] -0.90
#> [2,] -0.81
#>
#> $ab[[2]]
#> [,1] [,2]
#> [1,] -0.9 0.0
#> [2,] 0.0 -0.9
#>
#> $ab[[3]]
#> [,1] [,2] [,3]
#> [1,] -0.9 0.0 0
#> [2,] 0.0 -0.9 0
#>
#> $ab[[4]]
#> [,1] [,2] [,3] [,4]
#> [1,] -0.9 0.0 -2.614948e-16 3.228331e-16
#> [2,] 0.0 -0.9 0.000000e+00 0.000000e+00
#>
#>
alg1(t1[],c(5,5))
#> $orders
#> [1] 5 5
#>
#> $be
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 1 0.81 0.0000000 0 -4.343275e-16 -4.343275e-16
#> [2,] 1 0.90 0.6689647 0 0.000000e+00 7.861328e-16
#>
#> $fv
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 1 0.3439 0.3439000 0.3439000 0.3439000 0.3439000
#> [2,] 1 0.1900 0.1049724 0.1049724 0.1049724 0.1049724
#>
#> $bv
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 1 0.1900 0.19 0.19 0.19 0.19
#> [2,] 1 0.3439 0.19 0.19 0.19 0.19
#>
#> $af
#> $af[[1]]
#> [,1]
#> [1,] -0.81
#> [2,] -0.90
#>
#> $af[[2]]
#> [,1] [,2]
#> [1,] -0.8100000 0.0000000
#> [2,] -0.4972376 -0.4972376
#>
#> $af[[3]]
#> [,1] [,2] [,3]
#> [1,] -0.8100000 0.0000000 0
#> [2,] -0.4972376 -0.4972376 0
#>
#> $af[[4]]
#> [,1] [,2] [,3] [,4]
#> [1,] -0.8100000 0.0000000 -5.258951e-16 5.843279e-16
#> [2,] -0.4972376 -0.4972376 0.000000e+00 0.000000e+00
#>
#> $af[[5]]
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] -0.8100000 0.0000000 -1.05179e-15 5.843279e-16 5.843279e-16
#> [2,] -0.4972376 -0.4972376 0.00000e+00 5.258951e-16 -5.843279e-16
#>
#>
#> $ab
#> $ab[[1]]
#> [,1]
#> [1,] -0.90
#> [2,] -0.81
#>
#> $ab[[2]]
#> [,1] [,2]
#> [1,] -0.9 0.0
#> [2,] 0.0 -0.9
#>
#> $ab[[3]]
#> [,1] [,2] [,3]
#> [1,] -0.9 0.0 0
#> [2,] 0.0 -0.9 0
#>
#> $ab[[4]]
#> [,1] [,2] [,3] [,4]
#> [1,] -0.9 0.0 -2.614948e-16 3.228331e-16
#> [2,] 0.0 -0.9 0.000000e+00 0.000000e+00
#>
#> $ab[[5]]
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] -9.000000e-01 0.0 2.644003e-16 8.487282e-16 -1.057634e-15
#> [2,] 1.886404e-31 -0.9 0.000000e+00 -2.614948e-16 3.228331e-16
#>
#>