---
title: "Project 1_Time Series"
format:
html:
fig-format: png
code-fold: true
code-tools: true
embed-resources: true
mainfont: "Sukhumvit Set"
lang: th-TH
editor: visual
---
## Time Series Analysis: Visa (V) vs Mastercard (MA) - Q3 and Q4
นำเสนอโดย
1. นางสาวญาณิศา ฐิติกุล (รหัสนิสิต 6742059626)
2. นางสาวฐิตาภรณ์ ทิมศรี (รหัสนิสิต 6742026526)
3. นางสาวพิริยดา ลิ้มวัฒนายิ่งยง (รหัสนิสิต 6742077526)
------------------------------------------------------------------------
### ข้อ 1-2
#### Load packages and data
```{r}
library(quantmod)
library(fpp3) # tsibble, dplyr, tidyr, ggplot2, feasts, fable
```
```{r}
getSymbols("V", src = "yahoo", from = "2026-03-30", to = "2026-10-01")
getSymbols("MA", src = "yahoo", from = "2026-03-30", to = "2026-10-01")
range(zoo::index(V)); range(zoo::index(MA))
nrow(V); nrow(MA)
identical(zoo::index(V), zoo::index(MA))
```
สร้างตารางข้อมูลและ tsibble โดยใช้ลำดับวันทำการ (`day`) เป็น index เพราะวันที่จริงมีช่องว่างช่วงเสาร์-อาทิตย์
```{r}
prices <- inner_join(
data.frame(Date = zoo::index(V), V = as.numeric(Cl(V))),
data.frame(Date = zoo::index(MA), MA = as.numeric(Cl(MA))),
by = "Date"
) |> mutate(day = row_number())
write.csv(prices, "close_V_MA.csv", row.names = FALSE)
X <- prices |>
pivot_longer(c(V, MA), names_to = "Asset", values_to = "Close") |>
as_tsibble(index = day, key = Asset)
print(X)
```
------------------------------------------------------------------------
### ข้อ 3: Data Visualization
#### 3.1 กราฟแท่งเทียนและราคาปิดรายตัว
```{r}
chartSeries(V)
plot(V$V.Close, main = "V (Close Price in USD)")
chartSeries(MA)
plot(MA$MA.Close, main = "MA (Close Price in USD)")
```
#### 3.2 กราฟอนุกรมเวลาของทั้งสองสินทรัพย์
```{r}
autoplot(X, Close) + labs(y = "USD", x = "Trading day",title = "Closing price of V and MA")
plot.ts(prices$V, xlab = "Trading day", ylab = "USD", main = "V", las = 1)
plot.ts(prices$MA, xlab = "Trading day", ylab = "USD", main = "MA", las = 1)
```
#### 3.3 กราฟเปรียบเทียบแบบ normalized (วันแรก = 100) และ scatter plot
ราคาของสองตัวอยู่คนละระดับ จึงปรับให้เริ่มที่ 100 เพื่อเทียบการเคลื่อนไหว
```{r}
plot.ts(prices$V / prices$V[1] * 100, ylab = "Index (day 1 = 100)",
xlab = "Trading day", las = 1,
ylim = range(c(prices$V / prices$V[1], prices$MA / prices$MA[1]) * 100))
lines(prices$MA / prices$MA[1] * 100, col = "red")
legend("topleft", legend = c("V", "MA"), col = c("black", "red"), lty = 1)
library(GGally)
prices |> select(V, MA) |> ggpairs()
```
------------------------------------------------------------------------
### ข้อ 4: Trend, Seasonality, Cycle, Structural Change
#### 4.1 Trend: Centered Moving Average
ใช้ window 5 และ 21 วันทำการ (ประมาณ 1 สัปดาห์ และ 1 เดือน)
```{r}
for (nm in c("V", "MA")) {
s <- ts(prices[[nm]])
plot.ts(s, main = paste(nm, "- CMA"), xlab = "Trading day", ylab = "USD", las = 1)
lines(stats::filter(s, rep(1/5, 5), method = "convolution", sides = 2), lwd = 2, col = "red")
lines(stats::filter(s, rep(1/21, 21), method = "convolution", sides = 2), lwd = 2, col = "blue")
legend("topleft", col = c("red", "blue"), legend = c("5-MA", "21-MA"), lty = 1, lwd = 2)
}
```
#### 4.2 Trend เชิงปริมาณ: linear trend ด้วย TSLM
```{r}
X |> filter(Asset == "V") |> model(TSLM(log(Close) ~ trend())) |> report()
X |> filter(Asset == "MA") |> model(TSLM(log(Close) ~ trend())) |> report()
```
#### 4.3 Seasonality: STL decomposition
```{r}
X |> model(STL(log(Close) ~ season(period = 5, window = "periodic"))) |>
components() |> autoplot()
```
```{r}
y_v <- ts(log(prices$V), frequency = 5)
y_ma <- ts(log(prices$MA), frequency = 5)
plot(stl(y_v, s.window = "periodic"), main = "V: STL")
plot(stl(y_ma, s.window = "periodic"), main = "MA: STL")
```
#### 4.4 Classical Decomposition
```{r}
X |> model(classical_decomposition(Close ~ season(5), type = "multiplicative")) |>
components() |> autoplot()
```
#### 4.5 ตรวจ seasonality ด้วย ACF
```{r}
X |> ACF(log(Close), lag_max = 30) |> autoplot() + labs(title = "ACF of log(Close)")
X |> ACF(difference(log(Close)), lag_max = 30) |> autoplot() +
labs(title = "ACF of differenced log(Close)")
```
ทดสอบด้วย Ljung-Box ว่าการเปลี่ยนแปลงรายวัน (differenced log price) ยังมีสหสัมพันธ์เหลืออยู่หรือไม่
```{r}
for (nm in c("V", "MA")) {
d <- diff(log(prices[[nm]]))
cat("\n", nm, ": Ljung-Box (lag 10) ของ diff(log(Close))\n")
print(Box.test(d, lag = 10, type = "Ljung-Box"))
}
```
------------------------------------------------------------------------
### ข้อ 5:
#### Import Library
```{r}
library(quantmod)
```
```{r}
library(tseries)
```
### ข้อ 5.1 Descriptive analysis
Mastercard:
```{r}
# Define
x <- Cl(MA)
plot(x)
```
```{r}
plot.ts(MA$MA.Close, las=1, ylab="Mastercard Close")
```
```{r}
min(MA$MA.Close)
max(MA$MA.Close)
sd(MA$MA.Close)
```
```{r}
mean(MA$MA.Close) # sample mean
```
```{r}
var(MA$MA.Close) # sample var
```
```{r}
summary(MA$MA.Close)
```
Visa:
```{r}
plot.ts(V$V.Close, las=1, ylab="Visa Close")
```
```{r}
min(V$V.Close)
max(V$V.Close)
sd(V$V.Close)
```
```{r}
mean(V$V.Close) # sample mean
```
```{r}
var(V$V.Close) # sample var
```
```{r}
summary(V$V.Close)
```
### ข้อ 5.2 Stationary Test and White noise
```{r}
# Mastercard
# ADF Test
adf.test(MA$MA.Close)
```
p-value = 0.7289 \> 0.05 we fail to reject H0 -\> non-stationary
```{r}
# KPSS Test
kpss.test(MA$MA.Close)
```
p-value = 0.01 \< 0.05 we reject H0 -\> non-stationary
```{r}
# Visa
# ADF Test
adf.test(V$V.Close)
```
p-value = 0.6346 \> 0.05 we fail to reject H0 -\> non-stationary
```{r}
# KPSS Test
kpss.test(V$V.Close)
```
p-value = 0.01 \< 0.05 we reject the H0 -\> non-stationary
ดังนั้นทั้งในข้อมูลของ MA และ V ทั้งสอง test (ADF และ KPSS) ให้ผลลัพธ์ออกมาตรงกันว่าเป็น non-stationary -\> ต้องทำการ difference จนกว่าจะ stationary
#### ข้อ 5.3 1st order Differencing
```{r}
MA_dif1 <- na.omit(diff(log(Cl(MA))))
V_dif1 <- na.omit(diff(log(Cl(V))))
```
```{r}
# MA 1st order differencing
# ADF test
adf.test(MA_dif1)
# KPSS Test
kpss.test(MA_dif1)
```
in ADF Test: p-value \< 0.05 we reject H0 -\> stationary
in KPSS Test: p-value \> 0.05 we fail to reject H0 -\> stationary
```{r}
# V 1st order differencing
# ADF test
adf.test(V_dif1)
# KPSS Test
kpss.test(V_dif1)
```
in ADF Test: p-value \< 0.05 we reject H0 -\> stationary
in KPSS Test: p-value \> 0.05 we fail to reject H0 -\> stationary
ดังนั้นการทำ 1st order differencing ก็เพียงพอให้สามารถ stationary ได้ทั้งสำหรับ MA และ V
#### Descriptive stats ของ log return เพื่อดูค่าตอบแทน
```{r}
mean(MA_dif1) * 100
var(MA_dif1)
sd(MA_dif1) * 100
min(MA_dif1)
max(MA_dif1)
```
```{r}
mean(V_dif1) * 100
var(V_dif1)
sd(V_dif1) * 100
min(V_dif1)
max(V_dif1)
```
```{r}
# Sample autocovariance function
acf(MA_dif1, lag.max = 30, type = "covariance", plot = TRUE, las = 1)
```
```{r}
# Sample autocovariance at lag 0
acf(MA_dif1, type = "covariance", plot = FALSE)[0]
```
#### ข้อ 5.4 White noise
```{r}
# Sample autocorrelations
acf(MA_dif1, lag.max = 30, plot = FALSE)
```
```{r}
# MA - Box-Pierce Test (ไม่ใช้)
Box.test(MA_dif1, lag = 10, type = "Box-Pierce")
# MA - Ljung-Box test
Box.test(MA_dif1, lag = 10, type = "Ljung-Box")
```
from Box-Pierce, p-value \> 0.05 we fail to reject H0 -\> เป็น white noise
from Ljung-Box, p-value \> 0.05 we fail to reject H0 -\> เป็น white noise
```{r}
# V - Box-Pierce Test (ไม่ใช้)
Box.test(V_dif1, lag = 10, type = "Box-Pierce")
# V - Ljung-Box test
Box.test(V_dif1, lag = 10, type = "Ljung-Box")
```
from Box-Pierce, p-value \> 0.05 we fail to reject H0 -\> เป็น white noise
from Ljung-Box, p-value \> 0.05 we fail to reject H0 -\> เป็น white noise
------------------------------------------------------------------------
### ข้อ 6 ความสัมพันธ์เชิงเวลาระหว่างสินทรัพย์
#### ข้อ 6.1 sample correlation
```{r}
# sample correlation
asset <- na.omit(merge(MA_dif1, V_dif1))
cor(asset[,1], asset[, 2])
```
มี sample correlation อยู่ที่ 0.84 โดยประมาณ
#### ข้อ 6.2 CCF
```{r}
ccf(as.numeric(asset[,1]), as.numeric(asset[,2]), lag.max = 10,
main = "CCF of MA and V Log Returns")
```
```{r}
# ดูค่า lag0 จาก CCF
cc <- ccf(as.numeric(asset[,1]), as.numeric(asset[,2]),
lag.max = 10, plot = FALSE)
cc$acf[cc$lag == 0]
```
```{r}
1.96 / sqrt(nrow(asset)) # Bound
```
```{r}
# ดูค่า lag เป็นตาราง
data.frame(lag = cc$lag, ccf = round(cc$acf, 4))
```
เกินมาไม่เยอะมาก อาจเกิดจากข้อผิดพลาด
------------------------------------------------------------------------
### ข้อ 7 การวิเคราะห์แนวโน้มโดยใช้ตัวชี้วัดทางเทคนิค
#### SMA ของแต่ละสินทรัพย์
```{r}
#Plot SMA on V Close
plot.ts(V$V.Close, main = "SMA of V at Close Price", xlab = "Monthly", ylab = "Stock price",
xaxt = "n", las = 1)
dates <- seq(as.Date("2026-03-30"), as.Date("2026-10-01"), by = "day")
dates <- dates[!weekdays(dates) %in% c("Saturday", "Sunday")]
month_starts <- tapply(seq_along(dates), format(dates, "%Y-%m"), min)
month_starts <- month_starts[names(month_starts) != "2026-03"] # drop March
axis(1,
at = as.numeric(month_starts),
labels = format(dates[as.numeric(month_starts)], "%b"))
ma_5 <- stats::filter(V$V.Close, filter = rep(1/5,5), method = "convolution", sides = 1)
ma_21 <- stats::filter(V$V.Close, filter = rep(1/21,21), method = "convolution", sides = 1)
lines(ma_5, lwd = 2, col = "red") # 5-MA
lines(ma_21, lwd = 2, col = "blue") # 21-MA
legend("topleft", col = c("red", "blue"), legend = c("5-MA", "21-MA"),
lty = 1, lwd = 2)
```
```{r}
#Plot SMA on MA Close
plot.ts(MA$MA.Close, main = "SMA of MA at Close Price", xlab = "Monthly", ylab = "Stock price",
xaxt = "n", las = 1)
dates <- seq(as.Date("2026-03-30"), as.Date("2026-10-01"), by = "day")
dates <- dates[!weekdays(dates) %in% c("Saturday", "Sunday")]
month_starts <- tapply(seq_along(dates), format(dates, "%Y-%m"), min)
month_starts <- month_starts[names(month_starts) != "2026-03"] # drop March
axis(1,
at = as.numeric(month_starts),
labels = format(dates[as.numeric(month_starts)], "%b"))
ma_5 <- stats::filter(MA$MA.Close, filter = rep(1/5,5), method = "convolution", sides = 1)
ma_21 <- stats::filter(MA$MA.Close, filter = rep(1/21,21), method = "convolution", sides = 1)
lines(ma_5, lwd = 2, col = "red") # 5-MA
lines(ma_21, lwd = 2, col = "blue") # 21-MA
legend("topleft", col = c("red", "blue"), legend = c("5-MA", "21-MA"),
lty = 1, lwd = 2)
```
#### EMA ของแต่ละสินทรัพย์
```{r}
# EMA with fixed alpha (n=21) vs optimized alpha (min SSE)
n_fixed <- 21
alpha_fx <- 2 / (n_fixed + 1)
cat(sprintf("Fixed alpha for n=%d: %.4f\n\n", n_fixed, alpha_fx))
```
```{r}
# 1. EMA function (recursive definition)
# EMA_t = alpha * price_t + (1 - alpha) * EMA_{t-1}
# First value seeded with first price (common convention).
ema_recursive <- function(x, alpha) {
out <- numeric(length(x))
out[1] <- x[1]
for (i in 2:length(x)) {
out[i] <- alpha * x[i] + (1 - alpha) * out[i - 1]
}
return(out)
}
# 2. Optimize alpha by minimizing SSE (training set)
optimize_alpha <- function(train, alpha_grid = seq(0.01, 0.99, by = 0.005)) {
sse <- sapply(alpha_grid, function(a) {
ema <- ema_recursive(train, a)
sum((train - ema)^2, na.rm = TRUE)
})
best_alpha <- alpha_grid[which.min(sse)]
return(list(alpha = best_alpha, sse = min(sse), grid = alpha_grid, sse_all = sse))
}
# 3. Error metrics
metrics <- function(actual, pred) {
err <- actual - pred
c(
MAE = mean(abs(err), na.rm = TRUE),
RMSE = sqrt(mean(err^2, na.rm = TRUE)),
MAPE = mean(abs(err / actual), na.rm = TRUE) * 100
)
}
# 4. Run the analysis for one series
analyse_series <- function(price, label = "Series") {
price <- as.numeric(price)
n <- length(price)
cut <- floor(0.8 * n)
train <- price[1:cut]
test <- price[(cut + 1):n]
cat("=============================================\n")
cat(" Dataset:", label, " | n =", n, "| train =", cut, "| test =", n - cut, "\n")
cat("=============================================\n")
# -- (A) Fixed alpha = 2/(n+1) with n=21 --
ema_fx_full <- ema_recursive(price, alpha_fx)
pred_fx <- ema_fx_full[(cut + 1):n]
m_fx <- metrics(test, pred_fx)
# -- (B) Optimized alpha (min SSE on train) --
opt <- optimize_alpha(train)
ema_op_full <- ema_recursive(price, opt$alpha)
pred_op <- ema_op_full[(cut + 1):n]
m_op <- metrics(test, pred_op)
# -- (C) Baseline: simple trailing SMA(21) for reference --
sma21 <- stats::filter(price, rep(1 / 21, 21), sides = 1)
pred_sm <- sma21[(cut + 1):n]
m_sm <- metrics(test, pred_sm)
# -- Report --
cat(sprintf("Fixed alpha (a = %.4f, n = %d): MAE=%.4f RMSE=%.4f MAPE=%.2f%%\n",
alpha_fx, n_fixed, m_fx["MAE"], m_fx["RMSE"], m_fx["MAPE"]))
cat(sprintf("Optimized (a = %.4f, n_eq = %.2f): MAE=%.4f RMSE=%.4f MAPE=%.2f%%\n",
opt$alpha, 2 / opt$alpha - 1,
m_op["MAE"], m_op["RMSE"], m_op["MAPE"]))
cat(sprintf("SMA(21) reference: MAE=%.4f RMSE=%.4f MAPE=%.2f%%\n\n",
m_sm["MAE"], m_sm["RMSE"], m_sm["MAPE"]))
# -- Return everything for later use --
invisible(list(
fixed = list(alpha = alpha_fx, n = n_fixed, ema = ema_fx_full, metrics = m_fx),
optimized = list(alpha = opt$alpha, n_eq = 2 / opt$alpha - 1,
ema = ema_op_full, metrics = m_op,
grid = opt$grid, sse = opt$sse_all),
sma21 = list(ema = as.numeric(sma21), metrics = m_sm),
price = price,
train_cut = cut
))
}
# 5. Run on both datasets
res1 <- analyse_series(V$V.Close)
res2 <- analyse_series(MA$MA.Close)
```
#### Final Window Size Plot
#### SMA
```{r}
#SMA On V Close
plot.ts(V$V.Close, main = "SMA of V at Close Price", xlab = "Monthly", ylab = "Stock price",
xaxt = "n", las = 1)
dates <- seq(as.Date("2026-03-30"), as.Date("2026-10-01"), by = "day")
dates <- dates[!weekdays(dates) %in% c("Saturday", "Sunday")]
month_starts <- tapply(seq_along(dates), format(dates, "%Y-%m"), min)
month_starts <- month_starts[names(month_starts) != "2025-03"] # drop March
axis(1,
at = as.numeric(month_starts),
labels = format(dates[as.numeric(month_starts)], "%b"))
ma_21 <- stats::filter(V$V.Close, filter = rep(1/21,21), method = "convolution", sides = 1)
lines(ma_21, lwd = 2, col = "blue") # 21-MA
legend("topleft", col = "blue", legend = "21-MA", lwd = 2)
```
```{r}
#SMA On MA Close
plot.ts(MA$MA.Close, main = "SMA of MA at Close Price", xlab = "Monthly", ylab = "Stock price",
xaxt = "n", las = 1)
dates <- seq(as.Date("2025-03-30"), as.Date("2025-09-30"), by = "day")
dates <- dates[!weekdays(dates) %in% c("Saturday", "Sunday")]
month_starts <- tapply(seq_along(dates), format(dates, "%Y-%m"), min)
month_starts <- month_starts[names(month_starts) != "2025-03"] # drop March
axis(1,
at = as.numeric(month_starts),
labels = format(dates[as.numeric(month_starts)], "%b"))
ma_5 <- stats::filter(MA$MA.Close, filter = rep(1/5,5), method = "convolution", sides = 1)
ma_21 <- stats::filter(MA$MA.Close, filter = rep(1/21,21), method = "convolution", sides = 1)
lines(ma_21, lwd = 2, col = "blue") # 21-MA
legend("topleft", col = "blue", legend = "21-MA", lwd = 2)
```
#### EMA
```{r}
#EMA
price_V <- as.numeric(V$V.Close)
price_MA <- as.numeric(MA$MA.Close)
# Same length, same calendar (verify this!)
stopifnot(length(price_V) == length(price_MA))
n <- length(price_V)
# Numeric index (avoids the date-length problem)
idx <- 1:n
alpha_opt <- 0.99
ema_V_opt <- ema_recursive(price_V, alpha_opt)
ema_MA_opt <- ema_recursive(price_MA, alpha_opt)
```
```{r}
#Plot EMA comparing both stocks
plot(idx, ema_V_opt, type = "l", col = "blue", lwd = 2,
main = expression(paste("EMA (", alpha, " = 0.99) — Stock V vs Stock MA")),
xlab = "Trading day", ylab = "Price", las = 1,
ylim = range(c(ema_V_opt, ema_MA_opt), na.rm = TRUE))
lines(idx, ema_MA_opt, col = "red", lwd = 2)
legend("topleft", legend = c("Stock V", "Stock MA"),
col = c("blue","red"), lty = 1, lwd = 2, bty = "n")
norm_V <- ema_V_opt / ema_V_opt[1] * 100
norm_MA <- ema_MA_opt / ema_MA_opt[1] * 100
```
```{r}
#Normalized plot of EMA (starts at 100)
plot(idx, norm_V, type = "l", col = "blue", lwd = 2,
main = expression(paste("Normalized EMA (", alpha, " = 0.99)")),
xlab = "Trading day", ylab = "Index (Start = 100)", las = 1,
ylim = range(c(norm_V, norm_MA), na.rm = TRUE))
lines(idx, norm_MA, col = "red", lwd = 2)
abline(h = 100, lty = 3, col = "grey40")
legend("topleft", legend = c("Stock V", "Stock MA"),
col = c("blue","red"), lty = 1, lwd = 2, bty = "n")
```
