# load the ts dataset AirPassenger
data("AirPassengers")
# remove seasonality from a multiplicative model
AirPassengers_decomposed <- decompose(AirPassengers, type="multiplicative")
AirPassengers_seasonal_component <- AirPassengers_decomposed$seasonal
AirPassengers_seasonally_adjusted <- AirPassengers/AirPassengers_seasonal_component
par(mfrow=c(1,2))
plot.ts(AirPassengers, col = "blue", main = "Original series")
plot.ts(AirPassengers_seasonally_adjusted, col = "blue",
main = "Seasonally-adjusted series",
ylab = "Seasonally-adjusted values")
# remove the trend from a multiplicative model
AirPassengers_decomposed <- decompose(AirPassengers, type="multiplicative")
AirPassengers_trend_component <- AirPassengers_decomposed$trend
AirPassengers_detrended <- AirPassengers/AirPassengers_trend_component
par(mfrow=c(1,2))
plot.ts(AirPassengers, col = "blue", main = "Original series")
plot.ts(AirPassengers_detrended, col = "blue",
main = "Detrended series",
ylab = "Detrended values")
# create a simulate series
set.seed(1312)
toy_data <- arima.sim(n = 100, model = list(order = c(0,0,0)))
# add a deterministic trend to the series
toy_data_trend <- toy_data + 0.2*1:length(toy_data)
par(mfrow=c(1,3))
plot.ts(toy_data, main = "Original series")
plot.ts(toy_data_trend, main = "Series with Trend")
dummy_trend <- 1:length(toy_data_trend)
lm_toydata <- lm(toy_data_trend ~ dummy_trend)
plot.ts(lm_toydata$residuals, main = "Residuals (detrended)")
set.seed(111)
Random_Walk <- arima.sim(n = 500, model = list(order = c(0,1,0)))
Random_Walk_diff <- diff(Random_Walk)
par(mfrow=c(1,2))
plot.ts(Random_Walk,
main = "Random Walk",
col = "blue", ylab="")
plot.ts(Random_Walk_diff,
main = "Differenced Random Walk",
col = "blue", ylab="")
set.seed(123)
t <- 1:500
kappa <- 5 # costant term (or "drift")
delta <- 0.1
epsilon <- function(n){ # function for the error term
rnorm(n = 500, mean = 0, sd = 0.8)
}
y_I1 <- kappa + (delta * t) + arima.sim(n=499, list(order = c(0,1,0)),
rand.gen = epsilon)
y_I0 <- kappa + (delta * t) + arima.sim(n=500, list(order = c(1,0,0), ar = 0.8),
rand.gen = epsilon)
ts_y <- ts.intersect(y_I1, y_I0)
plot.ts(ts_y, plot.type = "single",
lty=c(1,3), col=c("red", "blue"),
main = "Stochastic w/ drift (red) Deterministic Trend (blue)",
ylab="")
# simulated data of x series correlated to y at lag 3 and 4
set.seed(999)
x_series <- arima.sim(n = 200, list(order = c(1,0,0), ar = 0.7, sd=1))
z <- ts.intersect(stats::lag(x_series, -3), stats::lag(x_series, -4))
y_series <- 15 + 0.8*z[,1] + 1.5*z[,2] + rnorm(196,0,1)
## Warning in `+.default`(15 + 0.8 * z[, 1] + 1.5 * z[, 2], rnorm(196, 0, 1)):
## longer object length is not a multiple of shorter object length
xy_series <- ts.intersect(y_series, z)
lm1 <- lm(xy_series[,1] ~ xy_series[,2] + xy_series[,3])
summary(lm1)
##
## Call:
## lm(formula = xy_series[, 1] ~ xy_series[, 2] + xy_series[, 3])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.57051 -0.73392 0.06238 0.74187 2.14794
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.96869 0.07238 206.80 <2e-16 ***
## xy_series[, 2] 0.85549 0.07680 11.14 <2e-16 ***
## xy_series[, 3] 1.42126 0.07678 18.51 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.002 on 196 degrees of freedom
## Multiple R-squared: 0.873, Adjusted R-squared: 0.8717
## F-statistic: 673.5 on 2 and 196 DF, p-value: < 2.2e-16
# install.package("forecast") # install the package if necessary
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
checkresiduals(lm1)
##
## Breusch-Godfrey test for serial correlation of order up to 10
##
## data: Residuals
## LM test = 8.689, df = 10, p-value = 0.5618
# another set of simulated data
# the x series is correlated at lag 3 and 4
set.seed(999)
x2_series <- arima.sim(n = 200, list(order = c(1,0,0), ar = 0.7, sd=1))
z2 <- ts.intersect(x2_series, stats::lag(x2_series, -3), stats::lag(x2_series, -4))
y2_series <- 15 + 0.8*z2[,2] + 1.5*z2[,3]
y2_errors <- arima.sim(n = 196, list(order = c(1,0,1), ar = 0.6, ma = 0.6), sd=1)
y2_series <- y2_series + y2_errors
# check the cross-correlations at lag 3 and 4
library(TSA)
## Registered S3 methods overwritten by 'TSA':
## method from
## fitted.Arima forecast
## plot.Arima forecast
##
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
##
## acf, arima
## The following object is masked from 'package:utils':
##
## tar
prew <- prewhiten(x2_series, y2_series)
prewhiten(x2_series, y2_series)
