/
residuals.R
193 lines (181 loc) · 5.99 KB
/
residuals.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
#' Residuals for various time series models
#'
#' Returns time series of residuals from a fitted model.
#'
#' Innovation residuals correspond to the white noise process that drives the
#' evolution of the time series model. Response residuals are the difference
#' between the observations and the fitted values (equivalent to \code{h}-step
#' forecasts). For functions with no \code{h} argument, \code{h=1}. For
#' homoscedastic models, the innovation residuals and the response residuals
#' for \code{h=1} are identical. Regression residuals are available for
#' regression models with ARIMA errors, and are equal to the original data
#' minus the effect of the regression variables. If there are no regression
#' variables, the errors will be identical to the original series (possibly
#' adjusted to have zero mean). \code{arima.errors} is a deprecated function
#' which is identical to \code{residuals.Arima(object, type="regression")}.
#' For \code{nnetar} objects, when \code{type="innovations"} and \code{lambda} is used, a
#' matrix of time-series consisting of the residuals from each of the fitted neural networks is returned.
#'
#' @param object An object containing a time series model of class \code{ar},
#' \code{Arima}, \code{bats}, \code{ets}, \code{arfima}, \code{nnetar} or
#' \code{stlm}.
#' If \code{object} is of class \code{forecast}, then the function will return
#' \code{object$residuals} if it exists, otherwise it returns the differences between
#' the observations and their fitted values.
#' @param type Type of residual.
#' @param h If \code{type='response'}, then the fitted values are computed for
#' \code{h}-step forecasts.
#' @param ... Other arguments not used.
#' @return A \code{ts} object.
#' @author Rob J Hyndman
#' @seealso \code{\link{fitted.Arima}}, \code{\link{checkresiduals}}.
#' @keywords ts
#'
#' @export
residuals.forecast <- function(object, type=c("innovation", "response"), ...) {
type <- match.arg(type)
if (type == "innovation") {
object$residuals
} else {
getResponse(object) - fitted(object)
}
}
#' @rdname residuals.forecast
#' @export
residuals.ar <- function(object, type=c("innovation", "response"), ...) {
type <- match.arg(type)
# innovation and response residuals are the same for AR models
object$resid
}
#' @rdname residuals.forecast
#'
#' @aliases residuals.forecast_ARIMA
#' @examples
#' fit <- Arima(lynx,order=c(4,0,0), lambda=0.5)
#'
#' plot(residuals(fit))
#' plot(residuals(fit, type='response'))
#' @export
residuals.Arima <- function(object, type=c("innovation", "response", "regression"), h=1, ...) {
type <- match.arg(type)
if (type == "innovation") {
object$residuals
} else if (type == "response") {
getResponse(object) - fitted(object, h = h)
} else {
x <- getResponse(object)
if (!is.null(object$lambda)) {
x <- BoxCox(x, object$lambda)
}
xreg <- getxreg(object)
# Remove intercept
if (is.element("intercept", names(object$coef))) {
xreg <- cbind(rep(1, length(x)), xreg)
}
# Return errors
if (is.null(xreg)) {
return(x)
} else {
norder <- sum(object$arma[1:4])
return(ts(
c(x - xreg %*% as.matrix(object$coef[(norder + 1):length(object$coef)])),
frequency = frequency(x), start = start(x)
))
}
}
}
#' @export
residuals.forecast_ARIMA <- residuals.Arima
#' @rdname residuals.forecast
#' @export
residuals.bats <- function(object, type=c("innovation", "response"), h=1, ...) {
type <- match.arg(type)
if (type == "innovation") {
object$errors
} else {
getResponse(object) - fitted(object, h = h)
}
}
#' @rdname residuals.forecast
#' @export
residuals.tbats <- function(object, type=c("innovation", "response"), h=1, ...) {
type <- match.arg(type)
if (type == "innovation") {
object$errors
} else {
getResponse(object) - fitted(object, h = h)
}
}
#' @rdname residuals.forecast
#' @export
residuals.ets <- function(object, type=c("innovation", "response"), h=1, ...) {
type <- match.arg(type)
if (type == "innovation") {
object$residuals
} else {
getResponse(object) - fitted(object, h = h)
}
}
#' @rdname residuals.forecast
#' @export
residuals.ARFIMA <- function(object, type=c("innovation", "response"), ...) {
type <- match.arg(type)
if (type == "innovation") {
if (!is.null(object$residuals)) { # Object produced by arfima()
return(object$residuals)
} else # Object produced by fracdiff()
{
if (is.element("x", names(object))) {
x <- object$x
} else {
x <- eval.parent(parse(text = as.character(object$call)[2]))
}
if (!is.null(object$lambda)) {
x <- BoxCox(x, object$lambda)
}
y <- fracdiff::diffseries(x - mean(x), d = object$d)
fit <- arima(y, order = c(length(object$ar), 0, length(object$ma)), include.mean = FALSE, fixed = c(object$ar, -object$ma))
return(residuals(fit, type = "innovation"))
}
}
else {
getResponse(object) - fitted(object)
}
}
#' @rdname residuals.forecast
#' @export
residuals.nnetar <- function(object, type=c("innovation", "response"), h=1, ...) {
type <- match.arg(type)
if (type == "innovation" && !is.null(object$lambda)) {
res <- matrix(unlist(lapply(object$model, residuals)), ncol = length(object$model))
if (!is.null(object$scalex$scale)) {
res <- res * object$scalex$scale
}
}
else {
res <- getResponse(object) - fitted(object, h = h)
}
tspx <- tsp(getResponse(object))
res <- ts(res, frequency = tspx[3L], end = tspx[2L])
return(res)
}
#' @rdname residuals.forecast
#' @export
residuals.stlm <- function(object, type=c("innovation", "response"), ...) {
type <- match.arg(type)
if (type == "innovation") {
object$residuals
} else {
getResponse(object) - fitted(object)
}
}
#' @rdname residuals.forecast
#' @export
residuals.tslm <- function(object, type=c("innovation", "response", "deviance"), ...) {
type <- match.arg(type)
if (type == "innovation" || type == "deviance") {
object$residuals
} else {
getResponse(object) - fitted(object)
}
}