Dynamic Linear Quantile Regression
dynrq(formula, tau = 0.5, data, subset, weights, na.action, method = "br", contrasts = NULL, start = NULL, end = NULL, ...)
formula |
a |
tau |
the quantile(s) to be estimated, may be vector valued, but all all values must be in (0,1). |
data |
an optional |
subset |
an optional vector specifying a subset of observations to be used in the fitting process. |
weights |
an optional vector of weights to be used
in the fitting process. If specified, weighted least squares is used
with weights |
na.action |
a function which indicates what should happen
when the data contain |
method |
the method to be used; for fitting, by default
|
contrasts |
an optional list. See the |
start |
start of the time period which should be used for fitting the model. |
end |
end of the time period which should be used for fitting the model. |
... |
additional arguments to be passed to the low level regression fitting functions. |
The interface and internals of dynrq are very similar to rq,
but currently dynrq offers two advantages over the direct use of
rq for time series applications of quantile regression:
extended formula processing, and preservation of time series attributes.
Both features have been shamelessly lifted from Achim Zeileis's
package dynlm.
For specifying the formula of the model to be fitted, there are several
functions available which allow for convenient specification
of dynamics (via d() and L()) or linear/cyclical patterns
(via trend(), season(), and harmon()).
These new formula functions require that their arguments are time
series objects (i.e., "ts" or "zoo").
Dynamic models: An example would be d(y) ~ L(y, 2), where
d(x, k) is diff(x, lag = k) and L(x, k) is
lag(x, lag = -k), note the difference in sign. The default
for k is in both cases 1. For L(), it
can also be vector-valued, e.g., y ~ L(y, 1:4).
Trends: y ~ trend(y) specifies a linear time trend where
(1:n)/freq is used by default as the covariate, n is the
number of observations and freq is the frequency of the series
(if any, otherwise freq = 1). Alternatively, trend(y, scale = FALSE)
would employ 1:n and time(y) would employ the original time index.
Seasonal/cyclical patterns: Seasonal patterns can be specified
via season(x, ref = NULL) and harmonic patterns via
harmon(x, order = 1). season(x, ref = NULL) creates a factor
with levels for each cycle of the season. Using
the ref argument, the reference level can be changed from the default
first level to any other. harmon(x, order = 1) creates a matrix of
regressors corresponding to cos(2 * o * pi * time(x)) and
sin(2 * o * pi * time(x)) where o is chosen from 1:order.
See below for examples.
Another aim of dynrq is to preserve
time series properties of the data. Explicit support is currently available
for "ts" and "zoo" series. Internally, the data is kept as a "zoo"
series and coerced back to "ts" if the original dependent variable was of
that class (and no internal NAs were created by the na.action).
###########################
## Dynamic Linear Quantile Regression Models ##
###########################
require(zoo)
## multiplicative median SARIMA(1,0,0)(1,0,0)_12 model fitted to UK seatbelt data
uk <- log10(UKDriverDeaths)
dfm <- dynrq(uk ~ L(uk, 1) + L(uk, 12))
dfm
dfm3 <- dynrq(uk ~ L(uk, 1) + L(uk, 12),tau = 1:3/4)
summary(dfm3)
## explicitly set start and end
dfm1 <- dynrq(uk ~ L(uk, 1) + L(uk, 12), start = c(1975, 1), end = c(1982, 12))
## remove lag 12
dfm0 <- update(dfm1, . ~ . - L(uk, 12))
tuk1 <- anova(dfm0, dfm1)
## add seasonal term
dfm1 <- dynrq(uk ~ 1, start = c(1975, 1), end = c(1982, 12))
dfm2 <- dynrq(uk ~ season(uk), start = c(1975, 1), end = c(1982, 12))
tuk2 <- anova(dfm1, dfm2)
## regression on multiple lags in a single L() call
dfm3 <- dynrq(uk ~ L(uk, c(1, 11, 12)), start = c(1975, 1), end = c(1982, 12))
anova(dfm1, dfm3)
###############################
## Time Series Decomposition ##
###############################
## airline data
## Not run:
ap <- log(AirPassengers)
fm <- dynrq(ap ~ trend(ap) + season(ap), tau = 1:4/5)
sfm <- summary(fm)
plot(sfm)
## End(Not run)
## Alternative time trend specifications:
## time(ap) 1949 + (0, 1, ..., 143)/12
## trend(ap) (1, 2, ..., 144)/12
## trend(ap, scale = FALSE) (1, 2, ..., 144)
###############################
## An Edgeworth (1886) Problem##
###############################
# DGP
fye <- function(n, m = 20){
a <- rep(0,n)
s <- sample(0:9, m, replace = TRUE)
a[1] <- sum(s)
for(i in 2:n){
s[sample(1:20,1)] <- sample(0:9,1)
a[i] <- sum(s)
}
zoo::zoo(a)
}
x <- fye(1000)
f <- dynrq(x ~ L(x,1))
plot(x,cex = .5, col = "red")
lines(fitted(f), col = "blue")Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.