Fit periodic time series models
fitPM.RdGeneric function with methods for fitting periodic time series models.
Arguments
- x
the time series.
- model
a periodic model, see Details.
- ...
further arguments to be passed on to individual methods.
Details
This is a generic function.
model provides the specification of the model. In particular,
the class of model determines what model is fitted. Specific
values of the parameters are generally ignored by non-iterative
methods but some methods can handle more detailed specifications, see
the individual methods.
Methods
signature(model = "ANY", x = "ANY")This is the default method. It simply exits with an error message stating that
fitPMdoes not have a method for the model specified bymodel.signature(model = "numeric", x = "ANY")Fits a PAR model to
x.modelshould be a vector of non-negative integers giving the PAR order. The length of this vector is taken to be the number of seasons.This is a convenience method. It constructs a PAR model and callls the method for
model = "PeriodicArModel".signature(model = "PeriodicArModel", x = "ANY")Fits a PAR model.
signature(model = "mcSpec", x = "ANY")Fits a periodic model according to the specification given by
model.Currently this method uses
mC.ssto set up the optimisation environment and then calls one of the optimisation functions in that environment as specified by argumentoptim.method, see below.Additional arguments may be specified to control the optimisation.
Argument
initcan be used to give initial values. It is passed on tomC.ss(and so has the format required by it).optim.methodis the name of an optimisation function in the environment returned bymC.ss. The default isoptim.method = "minim", which is based on the standard R functionoptim. Alternatives are "minimBB" or "minimBBLU". All this needs to be documented but seemC.ssandxx.ssfor details.Further arguments are passed on to the optimisation method. A typical argument supported by most optimisation functions is
control.signature(model = "PiPeriodicArModel", x = "ANY")Fits a periodically integrated PAR model using the parameters of
modelas initial values. Callspclspiarto do the actual work.signature(model = "SiPeriodicArModel", x = "ANY")Fits a seasonally integrated PAR model.
signature(model = "PeriodicArModel", x = "PeriodicMTS")signature(model = "PeriodicArModel", x = "PeriodicTS")
References
(todo: to be completed properly later)
Hipel KW, McLeod AI (1994). Time series modelling of water resources and environmental systems, Developments in water science; 45. London; Amsterdam: Elsevier.
Boshnakov GN, Iqelan BM (2009). “Generation of time series models with given spectral properties.” J. Time Series Anal., 30(3), 349--368. ISSN 0143-9782, doi: 10.1111/j.1467-9892.2009.00617.x .
Examples
## newm1 <- list(phi = matrix(1:12, nrow=4), p=rep(3,4), period=4, si2 = rep(1,4))
## new_pfm1 <- PeriodicFilterModel(newm1, intercept=0)
## generate some data;
set.seed(1234)
simts1 <- pcts(rnorm(1024), nseasons = 4)
fitPM(c(3,3,3,3), simts1)
#> An object of class FittedPeriodicArModel
#> Number of seasons: 4
#> Mean: -0.06076922 -0.002750496 -0.0114049 -0.05637432
#> SigmaSq: 1.060891 0.9714497 0.9903078 0.9074826
#>
#> Periodic order: AR(3,3,3,3)
#>
#> ar1 ar2 ar3
#> S1 0.11396728 -0.122469184 0.0075900305
#> S2 -0.00450794 0.083377197 0.0164112748
#> S3 -0.04491259 -0.005276924 0.0396226884
#> S4 0.03387766 -0.070385164 0.0001322766
#>
#> number of obs. for each season:
#> 256 256 256 256
fitPM(3, simts1)
#> An object of class FittedPeriodicArModel
#> Number of seasons: 4
#> Mean: -0.06076922 -0.002750496 -0.0114049 -0.05637432
#> SigmaSq: 1.060891 0.9714497 0.9903078 0.9074826
#>
#> Periodic order: AR(3,3,3,3)
#>
#> ar1 ar2 ar3
#> S1 0.11396728 -0.122469184 0.0075900305
#> S2 -0.00450794 0.083377197 0.0164112748
#> S3 -0.04491259 -0.005276924 0.0396226884
#> S4 0.03387766 -0.070385164 0.0001322766
#>
#> number of obs. for each season:
#> 256 256 256 256
## the fit on the underlying data is equivalent.
fitPM(c(3,3,3,3), as.numeric(simts1))
#> An object of class FittedPeriodicArModel
#> Number of seasons: 4
#> Mean: -0.06076922 -0.002750496 -0.0114049 -0.05637432
#> SigmaSq: 1.060891 0.9714497 0.9903078 0.9074826
#>
#> Periodic order: AR(3,3,3,3)
#>
#> ar1 ar2 ar3
#> S1 0.11396728 -0.122469184 0.0075900305
#> S2 -0.00450794 0.083377197 0.0164112748
#> S3 -0.04491259 -0.005276924 0.0396226884
#> S4 0.03387766 -0.070385164 0.0001322766
#>
#> number of obs. for each season:
#> 256 256 256 256
## equivalently, use a PAR(3,3,3,3) model for argument 'model'
## here the coefficients of pfm1 are ignored, since the estimation is linear.
pfm1 <- PeriodicArModel(matrix(1:12, nrow = 4), order = rep(3,4), sigma2 = 1)
pfm1
#> An object of class "PeriodicArModel"
#> Number of seasons: 4
#> Mean: 0 0 0 0
#> SigmaSq: 1 1 1 1
#>
#> Periodic order: AR(3,3,3,3)
#>
#> ar1 ar2 ar3
#> S1 1 5 9
#> S2 2 6 10
#> S3 3 7 11
#> S4 4 8 12
## these give same results as above
fitPM(pfm1, simts1)
#> An object of class FittedPeriodicArModel
#> Number of seasons: 4
#> Mean: -0.06076922 -0.002750496 -0.0114049 -0.05637432
#> SigmaSq: 1.060891 0.9714497 0.9903078 0.9074826
#>
#> Periodic order: AR(3,3,3,3)
#>
#> ar1 ar2 ar3
#> S1 0.11396728 -0.122469184 0.0075900305
#> S2 -0.00450794 0.083377197 0.0164112748
#> S3 -0.04491259 -0.005276924 0.0396226884
#> S4 0.03387766 -0.070385164 0.0001322766
#>
#> number of obs. for each season:
#> 256 256 256 256
fitPM(pfm1, as.numeric(simts1))
#> An object of class FittedPeriodicArModel
#> Number of seasons: 4
#> Mean: -0.06076922 -0.002750496 -0.0114049 -0.05637432
#> SigmaSq: 1.060891 0.9714497 0.9903078 0.9074826
#>
#> Periodic order: AR(3,3,3,3)
#>
#> ar1 ar2 ar3
#> S1 0.11396728 -0.122469184 0.0075900305
#> S2 -0.00450794 0.083377197 0.0164112748
#> S3 -0.04491259 -0.005276924 0.0396226884
#> S4 0.03387766 -0.070385164 0.0001322766
#>
#> number of obs. for each season:
#> 256 256 256 256
fitPM(c(1,1,1,1), simts1)
#> An object of class FittedPeriodicArModel
#> Number of seasons: 4
#> Mean: -0.06076922 -0.002750496 -0.0114049 -0.05637432
#> SigmaSq: 1.075896 0.978132 0.9917217 0.91232
#>
#> Periodic order: AR(1,1,1,1)
#>
#> ar1
#> S1 0.108486672
#> S2 0.001318166
#> S3 -0.041828251
#> S4 0.036776157
#>
#> number of obs. for each season:
#> 256 256 256 256
fitPM(c(3,2,2,1), simts1)
#> An object of class FittedPeriodicArModel
#> Number of seasons: 4
#> Mean: -0.06076922 -0.002750496 -0.0114049 -0.05637432
#> SigmaSq: 1.060947 0.9717132 0.9917187 0.91232
#>
#> Periodic order: AR(3,2,2,1)
#>
#> ar1 ar2 ar3
#> S1 0.113395504 -0.122778518 -0.0004315542
#> S2 -0.006365449 0.084234949 NA
#> S3 -0.041825811 -0.001666751 NA
#> S4 0.036776157 NA NA
#>
#> number of obs. for each season:
#> 256 256 256 256
fitPM(c(3,2,2,2), simts1)
#> An object of class FittedPeriodicArModel
#> Number of seasons: 4
#> Mean: -0.06076922 -0.002750496 -0.0114049 -0.05637432
#> SigmaSq: 1.060891 0.9717132 0.9917187 0.9074826
#>
#> Periodic order: AR(3,2,2,2)
#>
#> ar1 ar2 ar3
#> S1 0.113967280 -0.122469184 0.007590031
#> S2 -0.006365449 0.084234949 NA
#> S3 -0.041825811 -0.001666751 NA
#> S4 0.033877417 -0.070384980 NA
#>
#> number of obs. for each season:
#> 256 256 256 256
pdSafeParOrder(c(3,2,2,1))
#> [1] 3 2 2 2
pdSafeParOrder(rev(c(3,2,2,1)))
#> [1] 1 2 2 3
x <- arima.sim(list(ar = 0.9), n = 960)
pcx <- pcts(x, nseasons = 4)
mx <- matrix(x, nrow = 4)
##pc.acf(mx)
##pc.acf(mx, maxlag=10)
## TODO: avoid the warning when length ot the time series is not multiple
autocovariances(t(mx), maxlag = 6, nseasons = 4)
#> An object of class "Lagged2d"
#> Slot *data*:
#> Lag_0 Lag_1 Lag_2 Lag_3 Lag_4 Lag_5 Lag_6
#> 4.413696 3.698001 3.444325 3.083434 2.981296 2.740127 2.443547
#> 5.256342 4.480398 3.791234 3.376241 2.959491 2.945734 2.704898
#> 4.544553 4.226788 3.573730 2.787481 2.179323 1.972540 2.118425
#> 4.744290 4.149401 3.620436 3.039397 2.267672 1.595720 1.419777
#> 6.188442 4.913592 4.132992 3.413550 2.893992 1.944995 1.202820
#> 5.079685 5.095852 4.020852 3.462512 2.909262 2.517747 1.724930
#> 4.452634 4.097511 4.466851 3.581613 2.929054 2.568013 2.159905
#> 4.028403 3.766262 3.375314 3.895400 3.027806 2.462040 2.266185
#> 4.256483 3.789805 3.519559 3.181409 3.595804 2.723855 2.312791
#> 5.064014 4.082757 3.780803 3.857862 3.271179 3.849951 2.962659
#> 4.610889 4.241189 3.677709 3.378335 3.220566 2.679142 3.078473
#> 4.492845 4.083564 3.715123 3.132441 2.782178 2.780713 2.430421
#> 4.759835 4.221811 3.785254 3.334577 2.753388 2.379137 2.523684
#> 4.176337 3.914353 3.379820 3.318057 2.895464 2.089529 1.813878
#> 4.256341 3.795704 3.522276 3.120697 3.030260 2.548401 1.864357
#> 4.123065 3.730340 3.387307 3.057720 2.614952 2.406073 2.003338
autocovariances(t(mx))
#> An object of class "SampleAutocovariances"
#> , , Lag_0
#>
#>
#> 4.940244 4.425333 3.856613 3.312517
#> 4.425333 4.931782 4.153459 3.576585
#> 3.856613 4.153459 4.586589 4.037720
#> 3.312517 3.576585 4.037720 4.444731
#>
#> , , Lag_1
#>
#>
#> 3.053109 3.270405 3.765441 4.205794
#> 2.830618 3.009076 3.520014 3.771935
#> 2.262606 2.432067 2.839903 3.242108
#> 1.928944 2.007498 2.285101 2.672126
#>
#> , , Lag_2
#>
#>
#> 1.754047 1.767143 2.068452 2.409498
#> 1.430174 1.496301 1.828329 2.244051
#> 1.298815 1.204518 1.458011 1.818073
#> 1.238070 1.181317 1.385222 1.602964
#>
#> , , Lag_3
#>
#>
#> 1.3489948 1.2683819 1.4101905 1.573585
#> 0.9844681 0.8640092 1.0299877 1.269531
#> 0.9647423 0.9969533 0.9680465 1.158972
#> 1.1158544 1.1012862 1.0621654 1.144140
#>
#> , , Lag_4
#>
#>
#> 1.0681173 1.1595254 1.0799309 1.1913186
#> 0.8013338 0.9328923 0.7987529 0.8834352
#> 0.8382119 0.9644516 0.8820007 0.8776038
#> 0.9012098 1.0414786 1.0527977 0.9306610
#>
#> , , Lag_5
#>
#>
#> 0.8764845 1.0569533 1.0166630 0.9120474
#> 0.8499661 1.0786949 0.8387545 0.7325839
#> 0.6867834 0.9726713 0.8366954 0.8057214
#> 0.8406734 1.0936087 0.8225248 0.8450239
#>
#> , , Lag_6
#>
#>
#> 0.9444448 1.1035309 0.8113775 0.8291835
#> 1.1822682 1.2821839 1.0012574 0.8592283
#> 0.8407831 0.8361079 0.6282596 0.6109829
#> 0.8258035 0.7665833 0.6435341 0.6931841
#>
#> , , Lag_7
#>
#>
#> 0.8689502 0.7994249 0.7699847 0.7690557
#> 1.0024370 0.9676477 0.9242707 0.9797624
#> 0.6284261 0.5784530 0.5320890 0.5642380
#> 0.5484201 0.5434399 0.4715212 0.5167750
#>
#> , , Lag_8
#>
#>
#> 0.3468684 0.4347957 0.4606310 0.4804946
#> 0.6027137 0.5625705 0.6254156 0.6782344
#> 0.4822539 0.4593642 0.4413930 0.4271425
#> 0.4255004 0.4902419 0.4739370 0.3877285
#>
#> , , Lag_9
#>
#>
#> 0.2562200 0.3388748 0.2662562 0.2117790
#> 0.3822483 0.4504162 0.3251202 0.3949501
#> 0.1804663 0.3399479 0.2212054 0.3085179
#> 0.1899286 0.3206453 0.1264308 0.1618857
#>
#> , , Lag_10
#>
#>
#> 0.15180484 0.1875179 -0.03011971 -0.01144414
#> 0.05629723 0.1683812 -0.02502021 0.03307518
#> -0.08409547 0.1247324 -0.20909008 -0.16511840
#> -0.07978351 0.1421987 -0.13790670 -0.11767234
#>
#> , , Lag_11
#>
#>
#> -0.11875727 0.136673414 -0.12715441 -0.1241331
#> -0.13047205 0.071926993 -0.12785582 -0.1132037
#> -0.08908832 -0.009374151 -0.14856401 -0.1563577
#> -0.04172053 0.093328912 0.04497078 -0.1002592
#>
#> , , Lag_12
#>
#>
#> -0.0005362661 0.055885198 -0.06454129 -0.12617478
#> -0.0837065154 -0.030487010 -0.12492654 -0.05076665
#> -0.1603845066 -0.092526648 -0.18793024 -0.07192433
#> -0.0329581327 -0.002391903 -0.07344575 0.04800647
#>
#> , , Lag_13
#>
#>
#> 0.110795032 0.08497155 -0.10692939 0.039275700
#> 0.067135202 0.08472719 -0.08162884 -0.013675882
#> 0.003970013 0.07475564 -0.12919526 -0.199029101
#> 0.117818973 0.25836952 0.01765533 0.005315572
#>
#> , , Lag_14
#>
#>
#> 0.2569690 0.4167017 0.1750001 0.099989082
#> 0.2811720 0.3627738 0.2535473 0.142481793
#> 0.2898364 0.3371195 0.1645728 0.008941963
#> 0.2990899 0.3708633 0.1952323 0.025439383
#>
#> , , Lag_15
#>
#>
#> 0.3366081 0.3613924 0.3104424 0.1448003
#> 0.4606357 0.4791444 0.4525622 0.2297306
#> 0.3597465 0.3331262 0.3165522 0.1907237
#> 0.3781567 0.3046630 0.2260963 0.2302541
#>
#> , , Lag_16
#>
#>
#> 0.3386886 0.2301517 0.1967166 0.2457013
#> 0.4153102 0.2978708 0.3316045 0.3155847
#> 0.3095268 0.2480752 0.2709355 0.3258925
#> 0.2947141 0.2232702 0.2175240 0.3021335
#>
#> , , Lag_17
#>
#>
#> 0.3556765 0.2901412 0.2181435 0.3035475
#> 0.3473894 0.3594809 0.3021085 0.3914716
#> 0.3700382 0.4697207 0.3183025 0.4731177
#> 0.2458445 0.4474799 0.2438784 0.4462941
#>
#> Slot n:
#> [1] 240
#> Slot varnames:
#> [1] "Series 1" "Series 2" "Series 3" "Series 4"
#> Slot objectname: x
##It is an error to have more columns than rows.
## autocovariances(mx, maxlag = 6, nseasons = 4)
## autocovariances(mx)
num2pcpar(mx, c(1,1,1,1), period = 4)
#> $mean
#> [1] 0.2212892 0.2610529 0.2222639 0.1751942
#>
#> $coef
#> An object of class "slMatrix"
#> Slot "m":
#> lag
#> season 0 1
#> 1 1 -0.9462426
#> 2 1 -0.8957722
#> 3 1 -0.8421822
#> 4 1 -0.8803319
#>
#>
#> $scale
#> [1] 0.9800724 0.9837131 1.0433693 0.9435026
#>
num2pcpar(mx, c(3,3,3,3), period = 4)
#> $mean
#> [1] 0.2212892 0.2610529 0.2222639 0.1751942
#>
#> $coef
#> An object of class "slMatrix"
#> Slot "m":
#> lag
#> season 0 1 2 3
#> 1 1 -0.9962631 0.009726434 0.05118143
#> 2 1 -0.9003005 0.144883523 -0.15588601
#> 3 1 -0.7223822 -0.177308198 0.05138492
#> 4 1 -0.9444632 0.052370207 0.01986882
#>
#>
#> $scale
#> [1] 0.9768186 0.9723284 1.0339197 0.9404123
#>
sipfm1 <- new("SiPeriodicArModel", iorder = 1, siorder = 1, pcmodel = pfm1)
sipfm1
#> An object of class "SiPeriodicArModel"
#> iorder: 1, siorder: 1
#> PAR part of the model, An object of class "PeriodicArModel"
#> Number of seasons: 4
#> Mean: 0 0 0 0
#> SigmaSq: 1 1 1 1
#>
#> Periodic order: AR(3,3,3,3)
#>
#> ar1 ar2 ar3
#> S1 1 5 9
#> S2 2 6 10
#> S3 3 7 11
#> S4 4 8 12
fitPM(sipfm1, mx)
#> An object of class "SiPeriodicArModel"
#> iorder: 1, siorder: 1
#> PAR part of the model, An object of class FittedPeriodicArModel
#> Number of seasons: 4
#> Mean: 0 0 0 0
#> SigmaSq: 1.900178 2.014823 2.512701 1.919159
#>
#> Periodic order: AR(3,3,3,3)
#>
#> ar1 ar2 ar3
#> S1 0.03633172 0.02831852 0.02492281
#> S2 -0.02263934 -0.21176922 0.09351833
#> S3 -0.17347841 0.01862405 0.13451737
#> S4 0.02770379 0.11071311 0.07843985
#>
#> number of obs. for each season:
#> 238 238 238 238
pfm1
#> An object of class "PeriodicArModel"
#> Number of seasons: 4
#> Mean: 0 0 0 0
#> SigmaSq: 1 1 1 1
#>
#> Periodic order: AR(3,3,3,3)
#>
#> ar1 ar2 ar3
#> S1 1 5 9
#> S2 2 6 10
#> S3 3 7 11
#> S4 4 8 12
## experiments and testing
fit1 <- fitPM(c(3,3,3,3), simts1)
fit1_mf <- new("MultiFilter", coef = fit1@ar@coef)
vs <- mcompanion::mf_VSform(fit1_mf, form = "I")
tmp <- mcompanion::VAR2pcfilter(vs$Phi[ , -4],
Phi0inv = vs$Phi0inv, D = fit1@sigma2, what = "")
names(tmp) # "pcfilter" "var" "Uform"
#> [1] "pcfilter" "var" "Uform"
tmp$var
#> [1] 0.9074826 0.9903078 0.9714497 1.0608910
zapsmall(tmp$pcfilter)
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 0.11396728 -0.12246918 0.00759003 0 0 0
#> [2,] -0.00450794 0.08337720 0.01641127 0 0 0
#> [3,] -0.04491259 -0.00527692 0.03962269 0 0 0
#> [4,] 0.03387766 -0.07038516 0.00013228 0 0 0
fit1@ar@coef
#> lag
#> season 1 2 3
#> 1 0.11396728 -0.122469184 0.0075900305
#> 2 -0.00450794 0.083377197 0.0164112748
#> 3 -0.04491259 -0.005276924 0.0396226884
#> 4 0.03387766 -0.070385164 0.0001322766
all.equal(tmp$pcfilter[ , 1:3], fit1@ar@coef, check.attributes = FALSE) # TRUE
#> [1] TRUE
tmp$Uform
#> $Sigma
#> [1] 1.0608910 0.9714497 0.9903078 0.9074826
#>
#> $U0
#> [,1] [,2] [,3] [,4]
#> [1,] 1 -0.03387766 0.07038516 -0.0001322766
#> [2,] 0 1.00000000 0.04491259 0.0052769241
#> [3,] 0 0.00000000 1.00000000 0.0045079396
#> [4,] 0 0.00000000 0.00000000 1.0000000000
#>
#> $U
#> [,1] [,2] [,3]
#> [1,] 1.224810e-18 7.453890e-20 -2.117582e-22
#> [2,] 3.962269e-02 0.000000e+00 -6.776264e-21
#> [3,] 8.337720e-02 1.641127e-02 0.000000e+00
#> [4,] 1.139673e-01 -1.224692e-01 7.590031e-03
#>
#> $U0inv
#> [,1] [,2] [,3] [,4]
#> [1,] 1 0.03387766 -0.07190670 0.0002776578
#> [2,] 0 1.00000000 -0.04491259 -0.0050744608
#> [3,] 0 0.00000000 1.00000000 -0.0045079396
#> [4,] 0 0.00000000 0.00000000 1.0000000000
#>
#> $perm
#> [1] 4 3 2 1
#>
fit1@sigma2
#> [1] 1.0608910 0.9714497 0.9903078 0.9074826
## both give the matrix Sigma for the "I" form
identical(
vs$Phi0inv %*% diag(fit1@sigma2) %*% t(vs$Phi0inv)
,
tmp$Uform$U0inv %*% diag(tmp$Uform$Sigma) %*% t(tmp$Uform$U0inv)
) # TRUE
#> [1] TRUE
## no, this is a different matrix
var1_mat <- cbind(vs$Phi0, # identity matrix
- vs$Phi) # drop trailing zero columns?
var1_mat <- mcompanion::mCompanion(var1_mat)
var1_Sigma <- vs$Phi0inv %*% diag(fit1@sigma2) %*% t(vs$Phi0inv)
abs(eigen(diag(nrow(var1_mat)) - var1_mat)$values)
#> [1] 1.002444e+00 1.002444e+00 1.000000e+00 9.999852e-01 6.624797e-03
#> [6] 6.624797e-03 1.480313e-05 0.000000e+00