###############################
# Helper Functions
###############################
#' @importFrom scales percent
resid_hist = function(object){
ggplot(object, aes(x = residuals)) +
geom_histogram(aes(y=..density..), binwidth=.5) +
stat_density(geom="line", position="identity", aes(color = "Kernel")) +
stat_function(fun = dnorm, aes(color = "Normal")) +
#stat_function(fun = dt, args = list(df = df), aes(color = "Student's t")) +
theme_bw() + theme(legend.position =c(0.08, 0.85)) +
labs(
y = "Percent",
x = "Residuals",
colour = "Density"
) +
scale_colour_manual(values=c("grey", "blue", "orange")) +
scale_y_continuous(labels = scales::percent)
}
# Hack to use the gmwm gts plot object
resid_plot = function(object, ...) {
res = gts(object)
type = if(attr(object, "std")){ "Standardized Residuals" } else { "Residuals"}
plot(res) + ylab(type) + xlab("Observation Number")
}
#' Residuals of the Process
#'
#' Extracts the residuals of the process
#' @param x A \code{\link{arima}}, \code{\link{lm}}, or a data object.
#' @param std A \code{boolean} indicating whether the residuals should be
#' standardized.
#' @export
#' @rdname diag_resid
diag_resid = function(x, std = FALSE) {
UseMethod("diag_resid")
}
#' @export
#' @rdname diag_resid
diag_resid.Arima = function(x, std = FALSE) {
diag_resid.default(resid(x), std = std)
}
#' @export
#' @rdname diag_resid
diag_resid.lm = function(x, std = FALSE) {
diag_resid.default(resid(x), std = std)
}
#' @export
#' @rdname diag_resid
diag_resid.default = function(x, std = FALSE) {
if(is.ts(x)){
x = as.numeric(x)
}
if(std){
x = x/sd(x)
}
structure(data.frame(residuals = x, stringsAsFactors = FALSE),
std = std,
class = c("diag_resid","data.frame"))
}
#' Graph the Distribution of Residuals
#'
#' Makes a histogram of the distribution of residuals
#' @param x,object A \code{\link{diag_resid}} object.
#' @param type A \code{string} indicating either:
#' \code{"hist"} (residual histogram), \code{"resid"} (residual plot),
#' or \code{"both"}
#' @seealso \code{\link[gmwm]{gmwm}} and \code{\link[gmwm]{plot.gmwm}}
#' @return A \code{ggplot2} object.
#' @export
#' @rdname plot.diag_resid
plot.diag_resid = function(x, type = "hist", ...){
autoplot.diag_resid(object = x, type = type, ...)
}
#' @export
#' @rdname plot.diag_resid
autoplot.diag_resid = function(object, type = "hist", ...){
type = tolower(type)
if(type == "hist"){
resid_hist(object, ...)
} else if(type == "resid") {
resid_plot(object, ...)
} else{
grid.arrange(
resid_plot(object, ...),
resid_hist(object, ...),
nrow = 1
)
}
}
#' Wavelet Variance (WV) fit for White Noise
#'
#' Computes whether the residuals match with a
#' White Noise by using the Generalized Method of Wavelet Moments (GMWM)
#' @param x A \code{\link{arima}}, \code{\link{lm}}, or data set object.
#' @return A \code{diag_wv} object that inherits \code{gmwm} object information.
#' @export
#' @rdname diag_wv
diag_wv = function(x){
UseMethod("diag_wv")
}
#' @export
#' @rdname diag_wv
diag_wv.diag_resid = function(x){
diag_wv.default(resid(x))
}
#' @export
#' @rdname diag_wv
diag_wv.lm = function(x){
diag_wv.default(resid(x))
}
#' @export
#' @rdname diag_wv
diag_wv.Arima = function(x){
diag_wv.default(resid(x))
}
#' @export
#' @rdname diag_wv
diag_wv.default = function(x){
structure(gmwm(WN(),x), class = c("diag_wv","gmwm"))
}
#' Graph the Wavelet Variance (WV) Estimation of the White Noise
#'
#' Creates a plot containing the empirical WV and the implied WV for an
#' assumed white noise process.
#' @param x,object A \code{\link{diag_wv}} object.
#' @seealso \code{\link[gmwm]{gmwm}} and \code{\link[gmwm]{plot.gmwm}}
#' @return A \code{ggplot2} object.
#' @export
#' @rdname plot.diag_wv
plot.diag_wv = function(x, ...){
autoplot(x, ...)
}
#' @export
#' @rdname plot.diag_wv
autoplot.diag_wv = function(object, ...){
class(object) = "gmwm"
autoplot(object, title = "Haar WV") + theme(legend.position = "none")
}
#' QQ Normal plot functionality
#'
#' Provides the calculations behind a qq normal plot.
#' @inheritParams diag_resid
#' @return A \code{data.frame} with structure
#' \itemize{
#' \item{residuals}{Residuals}
#' \item{quantiles}{Quantiles}
#' \item{mirror}{Line y-values}
#' }
#' @export
#' @rdname diag_qq
diag_qq = function(x, std = FALSE){
UseMethod("diag_qq")
}
#' @export
#' @rdname diag_qq
#' @importFrom stats resid
diag_qq.Arima = function(x, std = FALSE){
diag_qq.default(resid(x), std = std)
}
#' @export
#' @rdname diag_qq
diag_qq.lm = function(x, std = FALSE){
diag_qq.default(resid(x), std = std)
}
#' @export
#' @rdname diag_qq
diag_qq.ts = function(x, std = FALSE){
diag_qq.default(data.matrix(x), std = std)
}
#' @export
#' @rdname diag_qq
diag_qq.default = function(x, std = FALSE){
if(is.ts(x)){
x = data.matrix(x)
}
if(std){
x = x/sd(x)
}
df_x = data.frame(residuals = sort(x))
n = nrow(df_x)
df_x$quantiles = qnorm((1:n)/(n+1))
m = (quantile(x, 0.75)-quantile(x,0.25))/
(qnorm(0.75)-qnorm(0.25))
b = quantile(x,0.25) - m*qnorm(0.25)
structure(df_x, m = m, b = b, std = std,
class = c("diag_qq","data.frame"))
}
#' QQ Normal Plot
#'
#' Generates a QQ Normal plot
#' @param x,object A \code{\link{diag_qq}} object
#' @param ... Additional parameters (not used)
#' @export
#' @rdname plot.diag_qq
plot.diag_qq = function(x, ...){
autoplot.diag_qq(object = x, ...)
}
#' @export
#' @rdname plot.diag_qq
autoplot.diag_qq = function(object, ...){
b = attr(object, "b")
m = attr(object, "m")
std = attr(object, "std")
class(object) = "data.frame"
ggplot(object, aes(x = quantiles)) +
geom_point(aes(y = residuals)) +
geom_abline(intercept = b,
slope = m,
color = "blue") +
labs(
x = "Theoretical Quantiles",
y = paste0(if(std){"Standardized "} else { "" }, "Sample Quantiles"),
title = "Normal Q-Q") +
theme_bw()
}
#' Bind Fitted and Residual Values
#'
#' Creates a \code{data.frame} containing fitted and residual values from a
#' model.
#' @param fitted A \code{numeric vector} containing fitted values from the model.
#' @param resid A \code{numeric vector} containing residuals.
#' @param std A \code{boolean} indicating whether to standardize the residuals
#' @return A \code{data.frame} with structure
#' \itemize{
#' \item{fitted}{Fitted Values}
#' \item{residuals}{Residuals}
#' }
#' @details
#' This function will be reworked at a later time to provide model support.
#' @export
#' @rdname diag_fitted
diag_fitted = function(fitted, resid, std = FALSE) {
if(std){
resid = resid/sd(resid)
}
structure(data.frame(fitted = as.numeric(fitted),
residuals = as.numeric(resid), stringsAsFactors = FALSE),
std = std,
class = c("diag_fitted","data.frame"))
}
#' Plot Fitted vs. Residual
#'
#' Generates a fitted vs. residual graph
#' @param x,object A \code{\link{diag_fitted}} object
#' @export
#' @rdname plot.diag_fitted
plot.diag_fitted = function(x, ...) {
autoplot.diag_fitted(object = x, ...)
}
#' @export
#' @rdname plot.diag_fitted
autoplot.diag_fitted = function(object, ...){
std = attr(object,'std')
restype = if(std){ "Standardized "} else { "" }
resabb = if(std) { "Std. "} else { "" }
ggplot(object, aes(x = fitted, y = residuals)) +
geom_point() +
geom_hline(yintercept = 0, color = "red", linetype="dashed") +
stat_smooth(method = "loess") +
labs(
x = "Fitted values",
y = paste0(restype, "Residuals"),
title = paste0(resabb, "Residuals vs Fitted Plot")
) + theme_bw()
}
#' Time Series Diagnostics
#'
#' Provides a list structure containing time series diagnostic information
#' @param model An \code{\link{arima}} object.
#' @param xt The data used to construct said model.
#' @param std A \code{boolean} indicating whether residuals should be standardized.
#' @param test A \code{string} indicating the type of portmanteau test to run.
#' @param robust A \code{boolean} indicating whether the ACF should be robust (\code{TRUE}) or
#' not (\code{FALSE}).
#' @return A \code{diag_ts} object
#' @export
diag_ts = function(model, xt, std = TRUE, test = "ljung-box", robust = FALSE){
res = resid(model)
fitted = as.numeric(xt - res) # fitted() no go for arima...
prt = if(tolower(test) == "ljung-box"){
diag_ljungbox(model, stdres = std, stop_lag = 30)
} else{
diag_boxpierce(model, stdres = std, stop_lag = 30)
}
if(std){
res2 = res/sd(res)
} else{
res2 = res
}
structure(list(
prt = prt,
res = diag_resid(model, std = std),
fit = diag_fitted(fitted, res, std = std),
acf = emp_acf(res2),
pacf = emp_pacf(res2),
wv = diag_wv(res2),
qq = diag_qq(model, std = std)), class = c("diag_ts","list"))
}
#' Graph Time Series Diagnostic
#'
#' Creates a time series diagnostic plot
#' @param x,object A \code{\link{diag_ts}} object.
#' @return A \code{grid.arrange} object.
#' @export
#' @rdname plot.diag_ts
plot.diag_ts = function(x, ...){
autoplot.diag_ts(object = x, ...)
}
#' @export
#' @rdname plot.diag_ts
autoplot.diag_ts = function(object, ...){
with(object,
grid.arrange(autoplot(res, type = "resid"), autoplot(fit), autoplot(res, type = "hist"), autoplot(qq),
autoplot(acf), autoplot(pacf),autoplot(wv), autoplot(prt), nrow = 2)
)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.