rm(list = ls(all = TRUE))
setwd("/Users/iyeonghan/skku/25-1/Time Series/Exercises")
## Import packages
library(itsmr)
library(TSlibrary)
##
## Attaching package: 'TSlibrary'
## The following objects are masked from 'package:itsmr':
##
## airpass, season, smooth.exp, smooth.ma, trend
library(tseries)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
##
## Attaching package: 'tseries'
## The following object is masked from 'package:itsmr':
##
## arma
library(urca)
mysterious <- read.table("~/skku/25-1/Time Series/Exercises/mysterious.txt", quote="\"", comment.char="")
df <- mysterious
remove(mysterious)
## Transform given data to ts type
data <- ts(df, start = c(2001, 1), end = c(2015, 12), frequency = 12)
## Time plot
plot.ts(data) ; title("Mysterious")
acf2(data, lag.max = 45); title("Correlogram of Data")
library(forecast)
##
## Attaching package: 'forecast'
## The following object is masked from 'package:itsmr':
##
## forecast
model <- tbats(data)
model$seasonal.periods
## NULL
##(b)
##Model1: Estimate trend by regression
x <- seq(from = 1, to = 180, by = 1)
x.sq <- x^2
x.cu <- x^3
out.lm <- lm(data ~ 1 + x + x.sq + x.cu)
summary(out.lm)
##
## Call:
## lm(formula = data ~ 1 + x + x.sq + x.cu)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.32439 -0.27213 -0.00626 0.32674 1.01982
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.398e-01 1.395e-01 -1.002 0.3176
## x 1.163e-02 6.654e-03 1.747 0.0823 .
## x.sq -1.720e-04 8.530e-05 -2.017 0.0452 *
## x.cu 6.298e-07 3.098e-07 2.033 0.0436 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.458 on 176 degrees of freedom
## Multiple R-squared: 0.03483, Adjusted R-squared: 0.01838
## F-statistic: 2.117 on 3 and 176 DF, p-value: 0.09975
lscurve <- ts(out.lm$fitted.values, start = c(2001,1), end = c(2015,12), frequency = 12)
plot.ts(data) ; title("Data and LS Curve")
lines(lscurve, col = "red")
model.1 <- data - lscurve
plot.ts(model.1); title("Result from Reg Based")
out.adf1 <- adf.test(model.1)
## Warning in adf.test(model.1): p-value smaller than printed p-value
print(out.adf1)
##
## Augmented Dickey-Fuller Test
##
## data: model.1
## Dickey-Fuller = -4.8119, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
out.adf1$p.value
## [1] 0.01
test(model.1)
## Null hypothesis: Residuals are iid noise.
## Test Distribution Statistic p-value
## Ljung-Box Q Q ~ chisq(20) 85.19 0 *
## McLeod-Li Q Q ~ chisq(20) 88.78 0 *
## Turning points T (T-118.7)/5.6 ~ N(0,1) 93 0 *
## Diff signs S (S-89.5)/3.9 ~ N(0,1) 90 0.8976
## Rank P (P-8055)/404.2 ~ N(0,1) 7926 0.7496
jarque.bera.test(model.1)
##
## Jarque Bera Test
##
## data: model.1
## X-squared = 0.61397, df = 2, p-value = 0.7357
## Model2: Estimate trend by smoothing
out.ex <- ts(smooth.exp(data, 0.4), start = c(2001,1), end = c(2015,12), frequency = 12)
plot.ts(data); lines(out.ex, col = "red")
title("Exponential Smoothing(0.4)")
model.2 <- data - out.ex
plot.ts(model.2); title("Result from Smoothing Based")
abline(h = 0, col = "red", lwd = 2)
out.adf2 <- ur.df(model.2, type="drift", selectlags="AIC")
summary(out.adf2)
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.75304 -0.15802 0.00131 0.16775 0.72661
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.003894 0.018188 0.214 0.831
## z.lag.1 -0.856880 0.081653 -10.494 < 2e-16 ***
## z.diff.lag 0.306746 0.071293 4.303 2.8e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2427 on 175 degrees of freedom
## Multiple R-squared: 0.3954, Adjusted R-squared: 0.3885
## F-statistic: 57.23 on 2 and 175 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -10.4942 55.079
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81
test(model.2)
## Null hypothesis: Residuals are iid noise.
## Test Distribution Statistic p-value
## Ljung-Box Q Q ~ chisq(20) 92.06 0 *
## McLeod-Li Q Q ~ chisq(20) 43.9 0.0016 *
## Turning points T (T-118.7)/5.6 ~ N(0,1) 93 0 *
## Diff signs S (S-89.5)/3.9 ~ N(0,1) 89 0.8976
## Rank P (P-8055)/404.2 ~ N(0,1) 8142 0.8296
jarque.bera.test(model.2)
##
## Jarque Bera Test
##
## data: model.2
## X-squared = 0.27533, df = 2, p-value = 0.8714
## Model3: Removing trend by Diff(1)
model.3 <- diff(data)
plot.ts(model.3); title("Differencing Based")
abline(h = 0, col = "red", lwd = 2)
out.adf3 <- ur.df(model.3, type="drift", selectlags="AIC")
summary(out.adf3)
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.33029 -0.28868 0.02541 0.27801 1.30953
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.005644 0.032766 0.172 0.8635
## z.lag.1 -1.210584 0.107641 -11.246 <2e-16 ***
## z.diff.lag 0.162083 0.074209 2.184 0.0303 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4359 on 174 degrees of freedom
## Multiple R-squared: 0.5334, Adjusted R-squared: 0.528
## F-statistic: 99.45 on 2 and 174 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -11.2465 63.2437
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81
test(model.3)
## Null hypothesis: Residuals are iid noise.
## Test Distribution Statistic p-value
## Ljung-Box Q Q ~ chisq(20) 57.01 0 *
## McLeod-Li Q Q ~ chisq(20) 32.42 0.039 *
## Turning points T (T-118)/5.6 ~ N(0,1) 114 0.476
## Diff signs S (S-89)/3.9 ~ N(0,1) 92 0.4386
## Rank P (P-7965.5)/400.8 ~ N(0,1) 8021 0.8899
jarque.bera.test(model.3)
##
## Jarque Bera Test
##
## data: model.3
## X-squared = 0.097821, df = 2, p-value = 0.9523
## Model4
out.ma12 <- smooth.ma(data, 12)
seasonal <- season(data - out.ma12, d = 12)
plot.ts(seasonal)
adj.data <- data - seasonal
plot.ts(adj.data)
model.4 <- diff(adj.data)
plot.ts(model.4)
acf(model.4, lag.max = 45)
out.adf4 <- adf.test(model.4)
## Warning in adf.test(model.4): p-value smaller than printed p-value
print(out.adf4)
##
## Augmented Dickey-Fuller Test
##
## data: model.4
## Dickey-Fuller = -8.9225, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
test(model.4)
## Null hypothesis: Residuals are iid noise.
## Test Distribution Statistic p-value
## Ljung-Box Q Q ~ chisq(20) 58.95 0 *
## McLeod-Li Q Q ~ chisq(20) 25.32 0.1895
## Turning points T (T-118)/5.6 ~ N(0,1) 114 0.476
## Diff signs S (S-89)/3.9 ~ N(0,1) 89 1
## Rank P (P-7965.5)/400.8 ~ N(0,1) 8036 0.8604
jarque.bera.test(model.4)
##
## Jarque Bera Test
##
## data: model.4
## X-squared = 0.48253, df = 2, p-value = 0.7856
library(knitr)
knitr::spin("2024mid.R", knit = FALSE)
## [1] "2024mid.Rmd"