Ekonometrik
Tugas 5
| Kontak | : \(\downarrow\) |
| yosia.yosia@student.matanauniversity.ac.id | |
| yyosia | |
| RPubs | https://rpubs.com/yosia/ |
Time Series
library(tidyquant)
library(ggplot2)
library(dplyr)
library(DT)
library(tsibble)
library(fpp3)
library(tidyverse)
library(fpp2)
library(tsbox)
library(tseries)
library(caTools)
library(forecast)Linear Model
multiple_stocks = read.csv("BBCA_JK.csv")
datatable(multiple_stocks)df = ts_ts(ts_long(multiple_stocks))
pricedata <- window(na.remove(df))
autoplot(pricedata[,c("High","Close")]) +
ylab("% change") + xlab("Year")
pricedata %>%
as.data.frame() %>%
ggplot(aes(x=Close, y=High)) +
ylab("High Price (quarterly % change)") +
xlab("Close Price (quarterly % change)") +
geom_point() +
geom_smooth(method="lm", se=FALSE)
tslm(High ~ Close, data=pricedata)##
## Call:
## tslm(formula = High ~ Close, data = pricedata)
##
## Coefficients:
## (Intercept) Close
## -11.850 1.049Fitted
fit.consMR <- tslm(
High ~ Open,
data=pricedata)
summary(fit.consMR)##
## Call:
## tslm(formula = High ~ Open, data = pricedata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -355.55 -101.34 -34.25 70.27 949.38
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11.234812 39.708220 -0.283 0.778
## Open 1.059689 0.009384 112.928 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 184.7 on 119 degrees of freedom
## Multiple R-squared: 0.9908, Adjusted R-squared: 0.9907
## F-statistic: 1.275e+04 on 1 and 119 DF, p-value: < 2.2e-16autoplot(pricedata[,'High'], series="Data") +
autolayer(fitted(fit.consMR), series="Fitted") +
xlab("Year") + ylab("") +
ggtitle("Percent change in BBCA High Price") +
guides(colour=guide_legend(title=" "))
cbind(Data = pricedata[,"High"],
Fitted = fitted(fit.consMR)) %>%
as.data.frame() %>%
ggplot(aes(x=Data, y=Fitted)) +
geom_point() +
ylab("Fitted (predicted values)") +
xlab("Data (actual values)") +
ggtitle("Percent change in BBCA High price") +
geom_abline(intercept=0, slope=1)
Forecasting with regression
Close2 <- window(pricedata[,'Close'], start=2018)
fit.Close <- tslm(Close2 ~ trend + season)
fcast <- forecast(fit.Close)
autoplot(fcast) +
ggtitle("Forecasts of Close Price using regression") +
xlab("Year") + ylab("megalitres")
summary(fit.Close)##
## Call:
## tslm(formula = Close2 ~ trend + season)
##
## Residuals:
## Min 1Q Median 3Q Max
## -703.96 -270.62 -12.71 271.00 786.46
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4984.160 227.421 21.916 < 2e-16 ***
## trend 48.994 4.551 10.765 8.36e-13 ***
## season2 -176.532 296.245 -0.596 0.5550
## season3 -535.526 295.930 -1.810 0.0787 .
## season4 -625.769 295.685 -2.116 0.0413 *
## season5 -629.763 295.510 -2.131 0.0400 *
## season6 -657.506 295.405 -2.226 0.0324 *
## season7 -445.250 295.370 -1.507 0.1404
## season8 -286.744 295.405 -0.971 0.3382
## season9 -476.987 295.510 -1.614 0.1152
## season10 -284.731 295.685 -0.963 0.3420
## season11 -162.474 295.930 -0.549 0.5864
## season12 34.782 296.245 0.117 0.9072
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 440.3 on 36 degrees of freedom
## Multiple R-squared: 0.8, Adjusted R-squared: 0.7334
## F-statistic: 12 on 12 and 36 DF, p-value: 3.086e-09Exponential Smoothing
Simple Exponential Smoothing
autoplot(pricedata[,("High")]) +
ylab("% change") + xlab("Year")
oildata <- window(pricedata[,"Close"], start= 2018)
# Estimate parameters
fc <- ses(oildata, h=5)
# Accuracy of one-step-ahead training errors
round(accuracy(fc),2)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 59.88 326.45 251.14 0.83 4.31 0.31 0.14autoplot(fc) +
autolayer(fitted(fc), series="Fitted") +
ylab("Close Price") + xlab("Year")
Trend Method
air <- window(pricedata[,"Open"], start= 2018)
fc <- holt(air, h=15)
fc2 <- holt(air, damped=TRUE, phi = 0.9, h=15)
autoplot(air) +
autolayer(fc, series="Holt's method", PI=FALSE) +
autolayer(fc2, series="Damped Holt's method", PI=FALSE) +
ggtitle("Forecasts from Holt's method") + xlab("Year") +
ylab("Open Price") +
guides(colour=guide_legend(title="Forecast"))
Estimation and model selection
aust <- window(pricedata[,"Close"], start=2018)
fit <- ets(aust)
summary(fit)## ETS(M,N,N)
##
## Call:
## ets(y = aust)
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 4530.8675
##
## sigma: 0.0572
##
## AIC AICc BIC
## 762.4352 762.9685 768.1106
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 63.15126 325.7903 248.4512 0.9066164 4.253759 0.3038405 0.1387965autoplot(fit)
Forecasting
fit %>% forecast(h=8) %>%
autoplot() +
ylab("International visitor night in Australia (millions)")
summary(fit)## ETS(M,N,N)
##
## Call:
## ets(y = aust)
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 4530.8675
##
## sigma: 0.0572
##
## AIC AICc BIC
## 762.4352 762.9685 768.1106
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 63.15126 325.7903 248.4512 0.9066164 4.253759 0.3038405 0.1387965Arima Models
indf_data <- getSymbols(Symbols = "BBCA.JK", src = "yahoo", from = Sys.Date() - 2016,
to = Sys.Date(), auto.assign = FALSE)
indf_data <- Cl(indf_data)growth rate of the rate
chart_Series(indf_data, col = "black")
add_SMA(n = 100, on = 1, col = "red")
add_SMA(n = 20, on = 1, col = "black")
add_RSI(n = 14, maType = "SMA")
add_BBands(n = 20, maType = "SMA", sd = 1, on = -1)
add_MACD(fast = 12, slow = 25, signal = 9, maType = "SMA", histogram = TRUE)
### Log tranformation stock data
indf_log <- log(indf_data)
head(indf_log, n = 10)## BBCA.JK.Close
## 2016-11-17 7.989560
## 2016-11-18 7.987864
## 2016-11-21 7.987864
## 2016-11-22 7.982758
## 2016-11-23 7.986165
## 2016-11-24 7.975908
## 2016-11-25 7.974189
## 2016-11-28 7.963808
## 2016-11-29 7.962067
## 2016-11-30 7.958577plot(indf_log, main = "log indf_data chart")
Karena data telah ditransformasikan secara log, kita dapat dengan jelas melihat bahwa rangkaian tersebut menunjukkan beberapa tren naik dan turun dalam interval waktu tertentu. Sepanjang tahun 2016 hingga pertengahan tahun 2017, saham menunjukkan tren kenaikan, sedangkan pada pertengahan tahun 2017 hingga awal tahun 2019 mengalami tren penurunan. Saham juga terdiri dari beberapa volatilitas dan ayunan. Inilah tanda-tanda pergerakan harga saham tidak stasioner.
Membedakan Data
### difference logged data
indf_diff <- diff(indf_log, lag = 1)
indf_diff <- na.locf(indf_diff, na.rm = TRUE,
fromLast = TRUE)
plot(indf_diff)
Augmented Dickey Fuller Test (Unit Root Test)
adf <- adf.test(indf_log, alternative = c("stationary", "explosive"),
k = 0)
adf##
## Augmented Dickey-Fuller Test
##
## data: indf_log
## Dickey-Fuller = -3.2859, Lag order = 0, p-value = 0.07311
## alternative hypothesis: stationaryadf_diff <- adf.test(indf_diff, alternative = c("stationary", "explosive"),
k = 0)
adf_diff##
## Augmented Dickey-Fuller Test
##
## data: indf_diff
## Dickey-Fuller = -39.805, Lag order = 0, p-value = 0.01
## alternative hypothesis: stationaryBuild the ARIMA model
train_data <- indf_diff[1:1384]
set.seed(123)
arima_model <- auto.arima(train_data, stationary = TRUE, ic = c("aicc", "aic", "bic"),
trace = TRUE)##
## Fitting models using approximations to speed things up...
##
## ARIMA(2,0,2) with non-zero mean : -7656.253
## ARIMA(0,0,0) with non-zero mean : -7648.948
## ARIMA(1,0,0) with non-zero mean : -7652.469
## ARIMA(0,0,1) with non-zero mean : -7653.953
## ARIMA(0,0,0) with zero mean : -7648.355
## ARIMA(1,0,2) with non-zero mean : -7655.207
## ARIMA(2,0,1) with non-zero mean : -7655.657
## ARIMA(3,0,2) with non-zero mean : -7654.646
## ARIMA(2,0,3) with non-zero mean : Inf
## ARIMA(1,0,1) with non-zero mean : -7651.478
## ARIMA(1,0,3) with non-zero mean : -7655.663
## ARIMA(3,0,1) with non-zero mean : -7654.096
## ARIMA(3,0,3) with non-zero mean : -7654.579
## ARIMA(2,0,2) with zero mean : -7655.526
##
## Now re-fitting the best model(s) without approximations...
##
## ARIMA(2,0,2) with non-zero mean : -7658.333
##
## Best model: ARIMA(2,0,2) with non-zero meansummary(arima_model)## Series: train_data
## ARIMA(2,0,2) with non-zero mean
##
## Coefficients:
## ar1 ar2 ma1 ma2 mean
## -1.1752 -0.5782 1.1208 0.4886 7e-04
## s.e. 0.2852 0.2543 0.3057 0.2760 4e-04
##
## sigma^2 = 0.0002302: log likelihood = 3835.2
## AIC=-7658.39 AICc=-7658.33 BIC=-7627
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -1.162948e-06 0.01514556 0.01006113 NaN Inf 0.6340565 -0.01013701checkresiduals(arima_model)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,0,2) with non-zero mean
## Q* = 9.8723, df = 5, p-value = 0.07894
##
## Model df: 5. Total lags used: 10Fitting ARIMA Model and Forecasting
arima <- arima(train_data, order = c(2, 0, 2))
summary(arima)##
## Call:
## arima(x = train_data, order = c(2, 0, 2))
##
## Coefficients:
## ar1 ar2 ma1 ma2 intercept
## -1.1752 -0.5782 1.1208 0.4886 7e-04
## s.e. 0.2852 0.2543 0.3057 0.2760 4e-04
##
## sigma^2 estimated as 0.0002294: log likelihood = 3835.2, aic = -7658.39
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -1.162948e-06 0.01514556 0.01006113 NaN Inf 0.6340565 -0.01013701forecast1 <- forecast(arima, h = 100)
plot(forecast1)
checkresiduals(forecast1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,0,2) with non-zero mean
## Q* = 9.8723, df = 5, p-value = 0.07894
##
## Model df: 5. Total lags used: 10arima <- arima(train_data, order = c(2,1, 2))
summary(arima)##
## Call:
## arima(x = train_data, order = c(2, 1, 2))
##
## Coefficients:
## ar1 ar2 ma1 ma2
## -0.7020 -0.1015 -0.3665 -0.6335
## s.e. 0.1389 0.0268 0.1377 0.1376
##
## sigma^2 estimated as 0.00023: log likelihood = 3826.79, aic = -7643.58
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.0002897375 0.01516129 0.0100553 NaN Inf 0.6336888 0.001936515forecast1 <- forecast(arima, h = 100)
plot(forecast1)
checkresiduals(forecast1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,2)
## Q* = 12.409, df = 6, p-value = 0.05344
##
## Model df: 4. Total lags used: 10