## Import libraries for analysis of time series
rm(list = ls(all = TRUE))
library(itsmr)
library(TSlibrary)
##
## Attaching package: 'TSlibrary'
## The following objects are masked from 'package:itsmr':
##
## airpass, season, smooth.exp, smooth.ma, trend
data = huron;
plot.ts(data);
title("Lake Huron Water level")
## Estimating global linear trend by regression
n = length(data);
x = seq(from = 1, to = n, by = 1);
out.lm1 = lm(data ~ 1 + x);
summary(out.lm1)
##
## Call:
## lm(formula = data ~ 1 + x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.50997 -0.72726 0.00083 0.74402 2.53565
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.202037 0.230111 44.335 < 2e-16 ***
## x -0.024201 0.004036 -5.996 3.55e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.13 on 96 degrees of freedom
## Multiple R-squared: 0.2725, Adjusted R-squared: 0.2649
## F-statistic: 35.95 on 1 and 96 DF, p-value: 3.545e-08
plot.ts(data);
title("Lake Huron Water level")
lines(out.lm1$fitted.values, col = "red")
plot.ts(data - out.lm1$fitted.values, main = "detrended", ylab = "detrended")
abline(h=0, col="red", lwd=2)
par(mfrow = c(2,2))
plot(out.lm1)
## Estimating global quadratic trend by regression
x_sq = seq(from = 1, to = n, by = 1)^2;
out.lm2 = lm(data ~ 1 + x + x_sq);
summary(out.lm2)
##
## Call:
## lm(formula = data ~ 1 + x + x_sq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6313 -0.6532 -0.1036 0.8176 2.5275
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.3165643 0.3169670 35.703 < 2e-16 ***
## x -0.0910728 0.0147787 -6.162 1.72e-08 ***
## x_sq 0.0006755 0.0001446 4.670 9.88e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.025 on 95 degrees of freedom
## Multiple R-squared: 0.4083, Adjusted R-squared: 0.3959
## F-statistic: 32.78 on 2 and 95 DF, p-value: 1.493e-11
par(mfrow = c(1,1))
plot.ts(data)
title("Lake Huron Water level")
lines(out.lm2$fitted.values, col = "red")
plot.ts(out.lm2$residuals, main = "Residuals", ylab = "Residuals")
abline(h=0, col="red", lwd=2)
par(mfrow = c(2,2))
plot(out.lm2)
## MA Filter(q=5)
par(mfrow = c(1,1))
out.ma = smooth.ma(data, 5)
plot.ts(data);
lines(out.ma, col = "red")
title("Detrend - MA")
plot.ts(data-out.ma);
title("Residuals - MA")
abline(h = 0, col = "red", lwd = 2)
## Exponential Smoothing
ex4 = smooth.exp(data, .4)
plot.ts(data, ylab = "level in feet")
title("Lake Huron Water Level")
lines(ex4, col = "red")
plot.ts(data - ex4);
title("Detrend - Exponential Smoothing")
abline(h = 0, col = "red", lwd = 2)
## Differencing (order = 1)
diff1 = diff(data)
diff2 = diff(diff1)
plot.ts(diff1)
title("Lag 1 Diff")
abline(h = 0, col = "red", lwd = 2)
## Differencing (order = 2)
plot.ts(diff2)
abline(h = 0, col = "red", lwd = 2)
title("Lag 2 Diff")
## Correlogram of Law Data
acf2(data, lag.max = 25)
title("Correlogram of Law Data")
## Correlogram of Residuals from linear trend fit
acf2(out.lm1$residuals, lag.max = 25)
title("Correlogram of Residuals from Linear Trend Fit")
## Test of Randomness for Detrended Data by First Order Regression
test(out.lm1$residuals)
## Null hypothesis: Residuals are iid noise.
## Test Distribution Statistic p-value
## Ljung-Box Q Q ~ chisq(20) 107.83 0 *
## McLeod-Li Q Q ~ chisq(20) 68.71 0 *
## Turning points T (T-64)/4.1 ~ N(0,1) 40 0 *
## Diff signs S (S-48.5)/2.9 ~ N(0,1) 50 0.6015
## Rank P (P-2376.5)/162.9 ~ N(0,1) 2344 0.8419
## Test of Randomness for Detrended Data by Second Order Regression
test(out.lm2$residuals)
## Null hypothesis: Residuals are iid noise.
## Test Distribution Statistic p-value
## Ljung-Box Q Q ~ chisq(20) 138.67 0 *
## McLeod-Li Q Q ~ chisq(20) 56.45 0 *
## Turning points T (T-64)/4.1 ~ N(0,1) 40 0 *
## Diff signs S (S-48.5)/2.9 ~ N(0,1) 50 0.6015
## Rank P (P-2376.5)/162.9 ~ N(0,1) 2406 0.8563
## Correlogram and Tests of Randomness of Residuals from MA Filter(q=5)
acf2(data - out.ma, lag.max = 25)
title("Correlogram of Detrended by MA(5)")
test(data - out.ma)
## Null hypothesis: Residuals are iid noise.
## Test Distribution Statistic p-value
## Ljung-Box Q Q ~ chisq(20) 130.2 0 *
## McLeod-Li Q Q ~ chisq(20) 24.88 0.2061
## Turning points T (T-64)/4.1 ~ N(0,1) 47 0 *
## Diff signs S (S-48.5)/2.9 ~ N(0,1) 46 0.3841
## Rank P (P-2376.5)/162.9 ~ N(0,1) 2355 0.895
## Correlogram and Tests of Randomness for Diff(1) and Diff(2)
acf2(diff1, lag.max = 25); title("Correlogram of Detrended by Diff(1)")
acf2(diff2, lag.max = 25); title("Correlogram of Detrended by Diff(2)")
test(diff1)
## Null hypothesis: Residuals are iid noise.
## Test Distribution Statistic p-value
## Ljung-Box Q Q ~ chisq(20) 25.24 0.1923
## McLeod-Li Q Q ~ chisq(20) 23.71 0.2554
## Turning points T (T-63.3)/4.1 ~ N(0,1) 63 0.9354
## Diff signs S (S-48)/2.9 ~ N(0,1) 48 1
## Rank P (P-2328)/160.4 ~ N(0,1) 2354 0.8713
test(diff2)
## Null hypothesis: Residuals are iid noise.
## Test Distribution Statistic p-value
## Ljung-Box Q Q ~ chisq(20) 23.74 0.2537
## McLeod-Li Q Q ~ chisq(20) 9.93 0.9694
## Turning points T (T-62.7)/4.1 ~ N(0,1) 63 0.9351
## Diff signs S (S-47.5)/2.8 ~ N(0,1) 49 0.5978
## Rank P (P-2280)/158 ~ N(0,1) 2349 0.6623
library(knitr)
knitr::spin("HW1.R", knit = FALSE)
## [1] "HW1.Rmd"