-
Notifications
You must be signed in to change notification settings - Fork 1
/
SMEDIAN.R
166 lines (135 loc) · 4.14 KB
/
SMEDIAN.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
#' @title Seasonal median model
#'
#' @description Train seasonal median model (SMEDIAN).
#'
#' @param .data Input data as tsibble.
#' @param specials Specials as list defined in \code{specials_smedian}.
#' @param ... Currently not in use.
#'
#' @return An object of class \code{SMEDIAN}.
train_smedian <- function(.data,
specials,
...){
if (length(measured_vars(.data)) > 1) {
abort("Only univariate responses are supported by SMEDIAN.")
}
y <- unclass(.data)[[measured_vars(.data)]]
if (all(is.na(y))) {
abort("All observations are missing, a model cannot be estimated without data.")
}
# Prepare period
lag <- specials$lag[[1]]
n <- length(y)
index <- rep(1:lag, times = ceiling(n / lag))[1:n]
smedian <- map_dbl(
.x = 1:lag,
.f = ~median(y[index == .x], na.rm = TRUE)
)
fitted <- rep(smedian, times = ceiling(n / lag))[1:n]
resid <- y - fitted
sigma <- sd(resid, na.rm = TRUE)
structure(
list(
fitted = fitted,
resid = resid,
smedian = smedian,
sigma = sigma,
last_period = last(index),
lag = lag),
class = "SMEDIAN"
)
}
globalVariables(c("self", "origin"))
specials_smedian <- new_specials(
lag = function(lag = NULL) {
lag <- get_frequencies(period = lag, data = self$data, .auto = "smallest")
if (lag == 1) {
abort("Non-seasonal model specification provided, use RW() or provide a different lag specification.")
}
if (!rlang::is_integerish(lag)) {
warn("Non-integer lag orders for random walk models are not supported. Rounding to the nearest integer.")
lag <- round(lag)
}
lag
},
.required_specials = c("lag")
)
#' @title Seasonal median model
#'
#' @description Automatic train a seasonal median model (SMEDIAN).
#'
#' @param formula Model specification (see "Specials" section, currently not in use ...).
#' @param ... Further arguments.
#'
#' @return smedian_model An object of class \code{SMEDIAN}.
#' @export
SMEDIAN <- function(formula, ...){
smedian_model <- new_model_class(
model = "SMEDIAN",
train = train_smedian,
specials = specials_smedian)
new_model_definition(
smedian_model,
!!enquo(formula),
...)
}
#' @title Forecast a trained seasonal median model
#'
#' @description Forecast a trained seasonal median model.
#'
#' @param object An object of class \code{SMEDIAN}.
#' @param new_data Forecast horizon (n-step ahead forecast)
#' @param specials Specials are currently not in use.
#' @param ... Additional arguments for forecast method.
#'
#' @return An object of class \code{fable}.
#' @export
forecast.SMEDIAN <- function(object,
new_data,
specials = NULL,
...){
# Extract model
smedian <- object$smedian
last_period <- object$last_period
n_ahead <- nrow(new_data)
index <- rep(1:length(smedian), times = ceiling((n_ahead + last_period) / length(smedian)))[(last_period + 1):(last_period + n_ahead)]
point <- as.numeric(smedian[index])
sd <- as.numeric(rep(object$sigma, times = n_ahead))
# Return forecasts
dist_normal(point, sd)
}
#' @title Extract fitted values from a trained seasonal median model
#'
#' @description Extract fitted values from a trained seasonal median model.
#'
#' @param object An object of class \code{SMEDIAN}.
#' @param ... Currently not in use.
#'
#' @return Fitted values as tsibble.
#' @export
fitted.SMEDIAN <- function(object, ...){
object[["fitted"]]
}
#' @title Extract residuals from a trained seasonal median model
#'
#' @description Extract residuals from a trained seasonal median model.
#'
#' @param object An object of class \code{SMEDIAN}.
#' @param ... Currently not in use.
#'
#' @return Fitted values as tsibble.
#' @export
residuals.SMEDIAN <- function(object, ...){
object[["resid"]]
}
#' @title Provide a succinct summary of a trained seasonal median model
#'
#' @description Provide a succinct summary of a trained seasonal median model.
#'
#' @param object An object of class \code{SMEDIAN}.
#'
#' @return Model summary as character value.
#' @export
model_sum.SMEDIAN <- function(object){
"SMEDIAN"
}