AR1 Transformer Function
AR1_transform.RdThe transformer function takes a bamlss.frame object and transforms
the response and the design matrices to account for lag 1 autocorrelation. The method
is also known as Prais-Winsten estimation.
trans_AR1(rho = 0.1)
AR1(rho = 0.1)Arguments
- rho
Specifies the correlation parameter at lag 1.
Value
A transformer function which can be used in the bamlss call.
References
Johnston, John (1972). Econometric Methods (2nd ed.). New York: McGraw-Hill. pp. 259--265.
See also
Examples
if (FALSE) ## Simulate AR1 data.
set.seed(111)
n <- 240
d <- data.frame("t" = 1:n)
## Nonlinear function.
f <- function(x) {
2 + sin(x / n * 2 * pi - pi)
}
## Correlated errors.
rho <- 0.8
e <- rnorm(n, sd = 0.1)
u <- c(e[1], rep(NA, n - 1))
for(i in 2:n){
u[i] <- rho * u[i - 1] + e[i]
}
## Response.
d$y <- f(d$t) + u
## Plot time-series data.
plot(d, type = "l")
## Estimate models without and with AR1 transformation.
b0 <- bamlss(y ~ s(t,k=20), data = d, criterion = "BIC")
#> BIC 136.4084 logPost -97.6646 logLik -44.2986 edf 8.7236 eps 0.3974 iteration 1
#> BIC -41.9715 logPost -8.4746 logLik 44.8913 edf 8.7236 eps 0.1812 iteration 2
#> BIC -152.399 logPost 46.7393 logLik 100.1053 edf 8.7236 eps 0.1157 iteration 3
#> BIC -189.094 logPost 75.9567 logLik 119.8168 edf 9.2214 eps 0.0649 iteration 4
#> BIC -221.762 logPost 91.1354 logLik 153.2188 edf 15.449 eps 0.0523 iteration 5
#> BIC -228.519 logPost 93.8649 logLik 162.6250 edf 17.649 eps 0.0209 iteration 6
#> BIC -228.766 logPost 94.0024 logLik 163.9364 edf 18.083 eps 0.0033 iteration 7
#> BIC -228.770 logPost 94.0090 logLik 164.1063 edf 18.144 eps 0.0002 iteration 8
#> BIC -228.770 logPost 94.0093 logLik 164.1187 edf 18.148 eps 0.0000 iteration 9
#> BIC -228.770 logPost 94.0093 logLik 164.1187 edf 18.148 eps 0.0000 iteration 9
#> elapsed time: 0.08sec
#> Starting the sampler...
#>
#> | | 0% 1.19sec
#> |* | 5% 2.94sec 0.15sec
#> |** | 10% 1.99sec 0.22sec
#> |*** | 15% 1.66sec 0.29sec
#> |**** | 20% 1.42sec 0.36sec
#> |***** | 25% 1.27sec 0.42sec
#> |****** | 30% 1.14sec 0.49sec
#> |******* | 35% 1.04sec 0.56sec
#> |******** | 40% 0.94sec 0.63sec
#> |********* | 45% 0.84sec 0.69sec
#> |********** | 50% 0.76sec 0.76sec
#> |*********** | 55% 0.68sec 0.83sec
#> |************ | 60% 0.59sec 0.89sec
#> |************* | 65% 0.52sec 0.96sec
#> |************** | 70% 0.44sec 1.02sec
#> |*************** | 75% 0.37sec 1.10sec
#> |**************** | 80% 0.29sec 1.17sec
#> |***************** | 85% 0.22sec 1.23sec
#> |****************** | 90% 0.14sec 1.30sec
#> |******************* | 95% 0.07sec 1.36sec
#> |********************| 100% 0.00sec 1.44sec
b1 <- bamlss(y ~ s(t,k=20), data = d, criterion = "BIC",
transform = AR1(rho = 0.8))
#> BIC 120.6978 logPost -100.564 logLik -19.7587 edf 14.812 eps 0.4482 iteration 1
#> BIC -183.664 logPost 54.4906 logLik 112.0547 edf 7.3797 eps 0.1809 iteration 2
#> BIC -316.134 logPost 120.7259 logLik 178.2900 edf 7.3797 eps 0.1060 iteration 3
#> BIC -372.549 logPost 148.9333 logLik 206.4974 edf 7.3797 eps 0.0648 iteration 4
#> BIC -381.197 logPost 153.2571 logLik 210.8212 edf 7.3797 eps 0.0266 iteration 5
#> BIC -380.659 logPost 169.9073 logLik 211.2022 edf 7.6169 eps 0.0059 iteration 6
#> BIC -380.659 logPost 169.9076 logLik 211.2025 edf 7.6169 eps 0.0002 iteration 7
#> BIC -380.659 logPost 169.9076 logLik 211.2025 edf 7.6169 eps 0.0000 iteration 8
#> BIC -380.659 logPost 169.9076 logLik 211.2025 edf 7.6169 eps 0.0000 iteration 8
#> elapsed time: 0.09sec
#> Starting the sampler...
#>
#> | | 0% 1.19sec
#> |* | 5% 1.20sec 0.06sec
#> |** | 10% 1.13sec 0.13sec
#> |*** | 15% 1.05sec 0.19sec
#> |**** | 20% 1.02sec 0.26sec
#> |***** | 25% 0.98sec 0.33sec
#> |****** | 30% 1.19sec 0.51sec
#> |******* | 35% 1.07sec 0.58sec
#> |******** | 40% 0.99sec 0.66sec
#> |********* | 45% 0.89sec 0.73sec
#> |********** | 50% 0.79sec 0.79sec
#> |*********** | 55% 0.71sec 0.87sec
#> |************ | 60% 0.62sec 0.93sec
#> |************* | 65% 0.54sec 1.01sec
#> |************** | 70% 0.46sec 1.07sec
#> |*************** | 75% 0.38sec 1.14sec
#> |**************** | 80% 0.30sec 1.21sec
#> |***************** | 85% 0.22sec 1.27sec
#> |****************** | 90% 0.15sec 1.34sec
#> |******************* | 95% 0.07sec 1.41sec
#> |********************| 100% 0.00sec 1.48sec
## Estimate full AR1 model.
b2 <- bamlss(y ~ s(t,k=20), data = d, criterion = "BIC",
family = "AR1")
#> BIC 128.8853 logPost -98.9894 logLik -37.7967 edf 9.7236 eps 0.5982 iteration 1
#> BIC -80.2651 logPost 5.5858 logLik 66.7784 edf 9.7236 eps 0.1685 iteration 2
#> BIC -246.753 logPost 88.8299 logLik 150.0225 edf 9.7236 eps 0.1025 iteration 3
#> BIC -339.747 logPost 135.3267 logLik 196.5194 edf 9.7236 eps 0.0643 iteration 4
#> BIC -363.204 logPost 147.0552 logLik 208.2479 edf 9.7236 eps 0.0326 iteration 5
#> BIC -348.719 logPost 157.7974 logLik 209.1956 edf 12.712 eps 0.0083 iteration 6
#> BIC -348.723 logPost 157.7992 logLik 209.1975 edf 12.712 eps 0.0006 iteration 7
#> BIC -348.723 logPost 157.7993 logLik 209.1975 edf 12.712 eps 0.0001 iteration 8
#> BIC -348.723 logPost 157.7993 logLik 209.1975 edf 12.712 eps 0.0000 iteration 9
#> BIC -348.723 logPost 157.7993 logLik 209.1975 edf 12.712 eps 0.0000 iteration 9
#> elapsed time: 0.15sec
#> Starting the sampler...
#>
#> | | 0% 1.67sec
#> |* | 5% 1.50sec 0.08sec
#> |** | 10% 1.49sec 0.17sec
#> |*** | 15% 1.43sec 0.25sec
#> |**** | 20% 1.38sec 0.35sec
#> |***** | 25% 1.33sec 0.44sec
#> |****** | 30% 1.26sec 0.54sec
#> |******* | 35% 1.19sec 0.64sec
#> |******** | 40% 1.10sec 0.73sec
#> |********* | 45% 1.01sec 0.83sec
#> |********** | 50% 0.93sec 0.93sec
#> |*********** | 55% 0.84sec 1.03sec
#> |************ | 60% 0.75sec 1.12sec
#> |************* | 65% 0.66sec 1.22sec
#> |************** | 70% 0.57sec 1.32sec
#> |*************** | 75% 0.47sec 1.42sec
#> |**************** | 80% 0.38sec 1.52sec
#> |***************** | 85% 0.28sec 1.61sec
#> |****************** | 90% 0.19sec 1.71sec
#> |******************* | 95% 0.10sec 1.81sec
#> |********************| 100% 0.00sec 1.90sec
rho <- predict(b2, model = "rho", type = "parameter")
print(range(rho))
#> [1] 0.7360218 0.7360218
## Estimated standard deviations.
sd0 <- predict(b0, model = "sigma", type = "parameter")
sd1 <- predict(b1, model = "sigma", type = "parameter")
sd2 <- predict(b2, model = "sigma", type = "parameter")
print(round(c(sd0[1], sd1[1], sd2[1]), 2))
#> [1] 0.13 0.00 0.10
## Plot fitted trends.
p0 <- predict(b0, model = "mu")
p1 <- predict(b1, model = "mu")
p2 <- predict(b2, model = "mu")
plot(d, type = "l")
lines(f(d$t) ~ d$t, col = 2, lwd = 2)
lines(p0 ~ d$t, col = 4, lwd = 2)
lines(p1 ~ d$t, col = 3, lwd = 3)
lines(p2 ~ d$t, col = 5, lwd = 3)
legend("topleft",
c("no trans", "with trans", "AR1 model", "truth"),
lwd = 2, col = c(4, 3, 5, 2), bty = "n")
