1.1. Một số package khai thác dữ liệu tài chính.

Trong R có một số package mà ta có thể dùng để khai thác dữ liệu tài chính như: tidyquant, tseries, quantmod. Các package này đều có một lượng đa dạng các hàm khai thác dữ liệu và mỗi một package đều hỗ trợ những kiểu dữ liệu khác nhau. Package tidyquant nằm trong hệ sinh thái tidyverse nên dữ liệu trả về của nó có dạng tibble object. Đây là một kiểu dữ liệu data frame nhẹ và không có chỉ số index đánh dấu mốc thời gian. Vì vậy khi hồi qui các model chuỗi thời gian ta cần phải convert sang định dạng chuỗi thời gian (zoo hoặc xts object). Tuy nhiên ưu điểm của tibble object đó là tích hợp các biến đổi của dplyr. Do đó các tính toán dữ liệu phức tạp được thực hiện dễ dàng hơn, đồng thời code được viết cũng gọn gàng hơn. Các package còn lại sẽ trả về dữ liệu dạng chuỗi thời gian và có thể sử dụng ngay vào các mô hình hồi qui timeseries.

• Lấy dữ liệu giá chứng khoán trên R sử dụng tidyquant, lưu ý định dạng trả về của tidyquant là tibble.

library(tidyquant)
GOOGLE1 <- tidyquant::tq_get('GOOGL', from = '2018-01-01', to = '2018-04-01')
head(GOOGLE1)

• Sử dụng tseries, dữ liệu trả về dạng chuỗi thời gian (zoo object).

library(tseries)
GOOGLE2 <- tseries::get.hist.quote(instrument = 'GOOGL', start = '2018-01-01', end = '2018-04-01')

## time series starts 2018-01-02
## time series ends   2018-03-29

head(GOOGLE2)

##               Open    High     Low   Close
## 2018-01-02 1053.02 1075.98 1053.02 1073.21
## 2018-01-03 1073.93 1096.10 1073.43 1091.52
## 2018-01-04 1097.09 1104.08 1094.26 1095.76
## 2018-01-05 1103.45 1113.58 1101.80 1110.29
## 2018-01-08 1111.00 1119.16 1110.00 1114.21
## 2018-01-09 1118.44 1118.44 1108.20 1112.79

• Sử dụng quanmod, dữ liệu trả về dạng chuỗi thời gian (xts object).

library(quantmod)
quantmod::getSymbols(Symbols = 'GOOGL', src = 'yahoo', from = '2018-01-01', to = '2018-04-01')

## [1] "GOOGL"

head(GOOGL)

##            GOOGL.Open GOOGL.High GOOGL.Low GOOGL.Close GOOGL.Volume
## 2018-01-02    1053.02    1075.98   1053.02     1073.21      1588300
## 2018-01-03    1073.93    1096.10   1073.43     1091.52      1565900
## 2018-01-04    1097.09    1104.08   1094.26     1095.76      1302600
## 2018-01-05    1103.45    1113.58   1101.80     1110.29      1512500
## 2018-01-08    1111.00    1119.16   1110.00     1114.21      1232200
## 2018-01-09    1118.44    1118.44   1108.20     1112.79      1340400
##            GOOGL.Adjusted
## 2018-01-02        1073.21
## 2018-01-03        1091.52
## 2018-01-04        1095.76
## 2018-01-05        1110.29
## 2018-01-08        1114.21
## 2018-01-09        1112.79

Điều đặc biệt của hàm getSymbols() trong quantmod đó là không cần phải lưu kết quả vào một tên bảng mà hàm số tự động export dữ liệu vào enviroment. Do đó ta không cần phải sử dụng phép gán GOOGL <- quantmod::getSymbols().

Muốn khai thác dữ liệu chứng khoán của Việt Nam có thể sử dụng VNDS package. Kết quả trả về định dạng tương tự như tidyquant.

library(VNDS)
VND <- VNDS::tq_get('VND', from = '2016-03-02', to = '2016-04-02')

## #VND from 2016-03-02 to 2016-04-02 already cloned

head(VND)

Vẽ biểu đồ chứng khoán: Bằng plot của package graphics, chỉ áp dụng đối với dữ liệu dạng xts hoặc zoo object.

GOOGLE <- GOOGLE2$Close
plot(GOOGLE, type="l")

Sử dụng ggplot2, áp dụng được cho cả xts object và data frame

library(ggplot2)
ggplot(GOOGLE1) +
  geom_line(aes(x = date, y = close))

1.2. Các định nghĩa trong time series.

Trung bình chuỗi: Xét chuỗi thười gian $x_t$. Trung bình của chuỗi được kí hiệu là $\mu_t$ và được tính như sau: $$\mu_t = \mathbf{E}(x_t)= \frac{\sum_{i = 1}^{N} x_i}{N}$$

Trong đó $\mathbf{E}$ là kí hiệu cho giá trị kì vọng của chuỗi.

Chuỗi trung bình trượt: Trung bình trượt là giá trị trung bình của một tập hợp các điểm thời gian liền kề. Mục đích của trung bình trượt là giảm nhiễu và tạo ra một chuỗi mượt hơn so với chuỗi gốc. Khi đó xu hướng biến đổi của chuỗi được thể hiện rõ hơn. Trong quá trình tính toán chuỗi trung bình trượt, ta sẽ sử dụng một cửa số với chiều dài $p$ và tiến hành trượt cửa sổ này dọc theo trục thời gian của chuỗi. Mỗi một đoạn trên chuỗi số mà cửa sổ trượt qua ta sẽ tính trung bình và thu được 1 điểm tương ứng thuộc chuỗi trung bình trượt. Độ dài của chuỗi chính là bậc của trung bình trượt. Trung bình trượt bậc p của chuỗi $x_t$ là chuỗi $v_t$ được xác định như sau: $$v_t = \frac{x_t+x_{t+1}+x_{t+2}+...+x_p}{p}$$

Bên dưới là ví dụ về đồ thị của một chuỗi trung bình trượt:

l = 200
x = sin(pi*seq(1, l, 1)/20) + rnorm(l, 0, 1)*0.2
plot(x, type = 'l')

x_mavg = c()
for (i in 1:(l-3)){
  x_mavg = c(x_mavg, mean(x[i:(i+3)]))
}
plot(x_mavg, type = 'l')

Ta có thể thấy trung bình trượt của chuỗi vẫn giữ được các đường nét chính của chuỗi gốc và đồng thời loại bỏ các nhiễu của chuỗi.

Hiệp phương sai: Với 2 chuỗi $x_t$ và $y_t$ bất kì làm thể làm để ta biết được 2 chuỗi này có mối quan hệ cùng chiều hay ngược chiều. Hiệp phương sai chuỗi sẽ có tác dụng đo lường mối quan hệ giữa 2 chuỗi số và được tính như sau:

$$\gamma_{x_t, y_t} = \mathcal{cov}(x_t, y_t) = \mathbf{E}[(x_t-\mathbf{E}(x_t))(y_t-\mathbf{E}(y_t))]$$

Khi $\gamma_{x_t, y_t} > 0$ hai chuỗi có mối quan hệ ngược chiều và trái lại. $\mathcal{cov}$ là viết tắt của covariance (hiệp phương sai).

Phương sai: Với một chuỗi $x_t$ cho trước ta luôn xác định được giá trị kì vọng của nó là $\mathbf{E}(x_t)$. Trung bình của tổng bình phương các độ lệch của chuỗi đối với kì vọng của chuỗi được gọi là phương sai:

$$\begin{eqnarray} \mathcal{var}(x_t) & = & \mathbf{E}[(x_t - \mathbf{E}(x_t)]^2 \\ & = & \mathbf{E}[x_t^2- 2x_t\mathbf{E}(x_t)+\mathbf{E}(x_t)^2] \\ & = & \mathbf{E}(x_t^2)-2\mathbf{E}[x_t\mathbf{E}(x_t)]+\mathbf{E}(x_t)^2 \\ & = & \mathbf{E}(x_t^2) - \mathbf{E}(x_t)^2 \end{eqnarray}$$

Phương sai là trường hợp đặc biệt của hiệp phương sai khi 2 chuỗi trùng nhau. $\mathcal{var}$ là kí hiệu viết tắt cho variance tức phương sai.

Độ lệch chuẩn: Ta nhận thấy rằng phương sai là trung bình của bình phương của độ lệch so với kì vọng. Do đó nếu ta lấy căn bậc hai của phương sai ta sẽ thu được đại lượng gọi là độ lệch chuẩn. Giá trị của độ lệch chuẩn được sử dụng để đánh giá mức độ biến động của dữ liệu.

$$\sigma(x_t) = \sqrt{\mathcal{var}(x_t)}$$

Hệ số tương quan: đo lường mối quan hệ giữa 2 chuỗi thời gian và có miền giá trị nằm trong khoảng [-1, 1]. Khi hệ số tương quan bằng 1, hai chuỗi số có mối quan hệ tuyến tính cùng chiều. Khi đó một biến có thể biểu diễn thông qua biến còn lại bằng 1 phương trình tuyến tính và hệ số góc của phương trình lớn hơn 0. Trái lại khi hệ số tương quan bằng -1 hai chuỗi có mối quan hệ tuyến tính là nghịch chiều và hệ số góc trong biểu diễn tuyến tính sẽ nhỏ hơn 0. Tại giá trị gần bằng 0, hệ số tương quan biểu thị cho 2 chuỗi độc lập với nhau. Sự biến đổi của chuỗi này không tác động lên chuỗi kia và ngược lại. Để tính được hệ số tương quan thì yêu cầu 2 chuỗi phải có số lượng phần tử bằng nhau. Hệ số tương quan của 2 chuỗi $x_t$ và $y_t$ được tính như sau: $$\rho_{x_t, y_t} = \frac{\mathbf{E}(x_t - \mathbf{E}(x_t))(y_t - \mathbf{E}(y_t))}{\sigma(x_t)\sigma(y_t)}$$ Tự tương quan: Tự tương quan là một khái niệm quan trọng trong chuỗi thời gian. Hầu hết các chuỗi thời gian sẽ có sự tương quan với giá trị trễ của nó và các giá trị càng gần nhau thì tương quan càng mạnh hoặc các giá trị cùng thuộc 1 chu kì của chuỗi thì sẽ có tương quan cao (chẳng hạn như cùng tháng trong chu kì năm hay cùng quí trong chu kì năm). Chính vì vậy hệ số này mới có tên là tự tương quan. Hệ số tự tương quan được viết tắt là ACF (auto correlation regression) và thường dùng để tìm ra độ trễ cần thiết để xây dựng các mô hình như ARIMA, GARCH, ARIMAX,… và kiểm tra yếu tố mùa vụ.

Hệ số tự tương quan bậc $s$ được xác định như sau:

$$\rho(s,t) = \frac{\gamma(s,t)}{\sqrt{\sigma_s\sigma_t}}$$

giá trị $\rho(s, t)$ đo lường khả năng dự báo của biến $x_t$ nếu chỉ sử dụng biến $x_s$. Trong trường hợp 2 đại lượng có tương quan hoàn hảo tức $\rho(s,t) = \pm 1$ ta có thể biểu diễn $x_t = \beta_0 + \beta_1 x_s$. Hệ số của $\beta_1$ sẽ ảnh hưởng lên chiều của hệ số tương quan. Theo đó $\rho(s, t) = 1$ khi $\beta_1 > 0$ và $\rho_(s, t) = -1$ khi $\beta_1 < 0$. Trong một số trường hợp ta cũng cần đo lường khả năng dự báo của biến $y_t$ được xác định bởi biến $x_s$. Khi đó hệ số tự tương quan sẽ không đáp ứng được do nó chỉ đo lường tường quan trên cùng chuỗi. Do đó chúng ta sẽ cần sử dụng đến hệ số tương quan chéo.

Hệ số tương quan chéo: (cross-correlation function) giữa 2 chuỗi $x_t$ và $y_t$ được xác định. $$\rho_{xy}(s, t) = \frac{\gamma_{xy}(s, t)}{\sqrt{\sigma_{x_s}\sigma_{y_t}}}$$ Nhiễu trắng (White noise): Một chuỗi $\epsilon_t$ được gọi là nhiễu trắng nếu nó thỏa mãn kì vọng bằng 0, các giá trị với các độ trễ khác nhau không có hiện tượng tự tương quan và phương sai sai số không đổi. Do kì vọng và phương sai không đổi nên chúng ta gọi phân phối của nhiễu trắng là phân phối xác định (identical distribution) và được kí hiệu là $\epsilon_t \sim WN(0, \sigma^2)$. Nhiễu trắng là một thành phần ngẫu nhiên thể hiện cho yếu tố không thể dự báo của model do nó không có tính qui luật. Trong R nhiễu trắng có thể tạo thành từ hàm arima.sim(). Bên dưới ta sẽ khởi tạo một chuỗi nhiễu trắng với độ dài 100.

wn = arima.sim(model = list(order = c(0, 0, 0)), n = 100)
plot(wn, type = 'l')

Ta nhận thấy nhiễu trắng có giá trị giao động ngẫu nhiên quanh giá trị 0 và đồng thời khoảng biến động không thay đổi theo thời gian. Đồ thị phân phối của chuỗi wn.

hist(wn, nclass = 10)

Đồ thị cho thấy wn có dạng phân phối chuẩn với kì vọng bằng 0. Kiểm tra tính tương quan của chuỗi wn với giá trị trễ bậc 1 của nó ta được:

library(ggplot2)
library(dplyr)
ggplot() + 
  geom_point(aes(x = wn, y = lag(c(wn), 1)))

Biểu đồ này cho thấy wn và lag(wn, 1) không có sự tương quan. Ngoài ra ta có thể sử dụng hàm acf() để biểu đồ hóa sự tự tương quan của các chuỗi với các độ trễ. Chúng ta gọi biểu đồ này là biểu đồ tự tương quan (autocorrelation graph):

acf(wn)

Các đường nét đứt màu xanh thể hiện cận trên và dưới của khoảng tin cậy 95% đối với kiểm định hệ số tương quan bằng 0. Trục tung là giá trị của hệ số tương quan trễ và trục hoành là các độ trễ tương ứng. Tại độ trễ bằng 0 hệ số tương quan đo lường chuỗi đó với chính nó nên giá trị bằng 1. Các độ trễ khác hệ số tương quan không vượt quá khoảng tin cậy 95% nên có thể thấy chuỗi nhiễu trắng không có hiện tượng tự tương quan.

Chuỗi dừng (stationary - timeseries): Là chuỗi không có xu hướng tăng hoặc giảm theo thời gian. Chúng ta biết rằng các ước lương và dự báo chuỗi thời gian đều dựa trên giả định chuỗi dừng. Tuy nhiên có rất nhiều chuỗi thời gian không đạt được tính dừng trên thực tế. Chẳng hạn như model với trend xác định.

$$y_t = c + \beta_1 t+\phi y_{t-1} + \theta \epsilon_{t-1} + \epsilon_{t}$$ Chúng ta cần chuyển qua chuỗi dừng bằng cách detrend thời gian để tạo thành một quá trình ngẫu nhiên có tính dừng:

$$y_t-\beta_1 t = c+\phi y_{t-1} + \theta \epsilon_{t-1} + \epsilon_{t}$$ Bước ngẫu nhiên: Bước ngẫu nhiên là chuỗi có tính chất: $$y_{t} = y_{t-1} + \epsilon_{t}$$ Với thành phần $\epsilon_t$ là nhiễu trắng, tức là: $\epsilon_t \sim WN(0, \sigma^2)$

Đồng tích hợp (intergration): Một chuỗi thời gian thông thường sẽ có 2 thành phần chính: Thành phần $TD_t$ được gọi là thành phần trend xác định (deterministic trend) đại diện cho yếu tố tương quan theo thời gian của chuỗi khiến chuỗi tăng hoặc giảm dần qua thời gian và thành phần còn lại có tính chất tự tương quan. Chẳng hạn chuỗi $y_t$ được biểu diễn theo 2 thành phần như bên dưới: $$\begin{eqnarray}y_t & = & TD_t+ z_t \\ TD_t & = & \alpha + \delta t \\ z_t & = & \phi z_{t-1} + \epsilon_t, \epsilon_t \sim WN(0, \sigma^2) \end{eqnarray}$$ $z_t$ là quá trình tự hồi qui bậc 1 hay còn gọi là AR(1).

Khi đó có 2 khả năng xảy ra ảnh hưởng tới hình dạng của $y_t$: • Nếu $|\phi| < 1$: $y_t$ là quá trình $\mathbf{I}(0)$ (đồng tích hợp bậc 0) đối với thành phần $TD_t$. • Nếu $\phi = 1$ thì $z_t = z_{t-1} + \epsilon_t = z_0 + \sum_{i=1}^{t}\epsilon_{i}$ là một quá trình ngẫu nhiên và $y_t$ là quá trình $\mathbf{I}(1)$ (đồng tích hợp bậc 1) đối với thành phần $TD_t$ với thành phần đuôi (drift).

Giả sử $\alpha = 5, \delta = 0.2$. Bên dưới ta sẽ mô phỏng hình dạng của 2 chuỗi cho 2 trường hợp:

• Trường hợp $\phi = 0.8$ ta có:

n = 250
t = seq(1, n, 1)
al = 5
del = 0.2
phi = 0.8
epl = rnorm(n, 0, 1)

TD = al + del*t
z = c(1) 
for (i in 1:(n-1)){
  z_t = z[length(z)]*phi + epl[i]
  z = c(z, z_t)
}

y1 = TD + z

• Trường hợp $\phi = 1$:

n = 250
t = seq(1, n, 1)
al = 5
del = 0.2
phi = 1
epl = rnorm(n, 0, 1)

TD = al + del*t
z = c(1) 
for (i in 1:(n-1)){
  z_t = z[length(z)]*phi + epl[i]
  z = c(z, z_t)
}

y2 = TD + z

plot(y1, type = 'l') 
lines(y2, lty = 2)

Ta nhận thấy $y_t$ có xu hướng tăng và xu hướng này bám sát hình thái tăng ban đầu của yếu tố trend $\mathbf{TD}_{t}$. Tuy nhiên trường hợp $\phi = 1$ thì $y_t$ vẫn có xu hướng tăng nhưng cao hơn so với yếu tố trend. Như vậy đồng tích hợp là quá trình cùng tăng hoặc cùng giảm của các chuỗi thời gian theo cũng một trend.

1.3. Kiểm tra tính dừng

Giả sử ta có một quá trình tự hồi qui AR(1) đối với chuỗi $y_t$ được xác định như sau: $$y_t = \phi y_{t-1} + \epsilon_t, \epsilon_t \sim WN(0, \sigma^2)$$ Khi khai triển $y_t$ một cách liên tục theo giá trị trễ ta có: $$\begin{eqnarray} y_t& = & \phi y_{t-1} + \epsilon_t \\ & = & \phi(\phi y_{t-2} + \epsilon_{t-1}) + \epsilon_t \\ & = &\dots \\ & = & \phi(\phi( \dots (\phi y_{0} + \epsilon_{0})+ \dots) + \epsilon_{t-1}) + \epsilon_t \\ & = & \phi^{t}y_0 + \phi^{t-1}\epsilon_{1} + \dots + \phi \epsilon_{t-1} + \epsilon_t \end{eqnarray}$$

Do đó $$\mathbf{E}(y_t) = \phi^{t}\mathbf{E}(y_0)+\mathbf{E}(\phi^{t-1}\epsilon_{1} + \dots + \phi \epsilon_{t-1} + \epsilon_t) = \phi^{t}\mathbf{E}(y_0)$$ Dấu bằng thứ 2 xảy ra là do $\mathbf{E}(\epsilon_t) = 0, \forall t \in \mathbb{N}$ Mặt khác $$\begin{equation} \begin{cases} \phi > 1: \lim\limits_{t \rightarrow \infty} \phi^{t} \mathbf{E}(y_0) = \infty \\ \phi = 1: \lim\limits_{t \rightarrow \infty} \phi^{t} \mathbf{E}(y_0) = \mathbf{E}(y_0) \\ \phi < 1: \lim\limits_{t \rightarrow \infty} \phi^{t} \mathbf{E}(y_0) = 0 \end{cases} \end{equation}$$ Do đó chuỗi $y_t$ sẽ dừng khi $\phi \leq 1$. Như vậy kiểm định nghiệm đơn vị sẽ dựa trên cặp giả thuyết: $$\begin{equation} \begin{cases} H_0: \phi = 1, \Rightarrow y_t \sim \mathbf{I}(1) \\ H1: |\phi| < 1, \Rightarrow y_t \sim \mathbf{I}(0) \end{cases} \end{equation}$$ Chúng ta gọi đó là kiểm định nghiệm đơn vị (unit root test) bởi vì phương trình đặc trưng của nó $\theta(y) = 1-\theta L = 0$ có nghiệm đơn vị. Giá trị ngưỡng kiểm định:

$$DF = \frac{\hat\phi - 1}{SE(\hat\phi)}$$

Chúng ta sẽ so sánh giá trị ngưỡng kiểm định này với giá trị tới hạn của phân phối Dickey - Fuller để đưa ra kết luận về chấp nhận hoặc bác bỏ giả thuyết $H_0$. Kiểm định tính dừng của chuỗi, sử dụng hàm adf.test() của package tseries. Định dạng đầu vào phải có dạng vector:

library(tseries)
plot(GOOGLE)

adf.test(GOOGLE)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  GOOGLE
## Dickey-Fuller = -3.295, Lag order = 3, p-value = 0.08085
## alternative hypothesis: stationary

Trong các kiểm định thống kê của R, hầu hết các kết quả trả về đều có dạng xác xuất. Việc này sẽ thuận tiện hơn cho người dùng khi đưa ra kết luận vì ngưỡng tới hạn để bác bỏ giả thuyết luôn là 0.05. Do đó thay vì so sánh cặp giá trị kiểm định với giá trị tới hạn ta sẽ so sánh giá trị xác xuất kiểm định p-value với ngưỡng 0.05. Khi đã đọc nhiều kiểm định sẽ hình thành phản xạ rất nhanh khi so sánh với 0.05. Giả thuyết $H_0$ tức chuỗi không dừng hoặc đồng tích hợp bậc 1 được chấp nhận khi p-value > 0.05 và chuỗi dừng khi p-value < 0.05.

Trong tính huống này ta kết luận chuỗi không có tính dừng vì kiểm định có p-value > 0.05. Lấy sai phân trong tình huống chuỗi không dừng.

GOOGLE.DIFF <- diff(GOOGLE)

Kiểm tra tính dừng

plot(GOOGLE.DIFF, type = 'l')

adf.test(GOOGLE.DIFF)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  GOOGLE.DIFF
## Dickey-Fuller = -2.719, Lag order = 3, p-value = 0.2842
## alternative hypothesis: stationary

Bên cạnh sử dụng sai phân chuỗi ta cũng có thể dùng sai phân của logarit để tạo chuỗi dừng. Sai phân của logarit được lựa chọn là do nó thể hiện phần trăm tăng trưởng của chuỗi qua thời gian.

LOG.GOOGLE.DIFF <- diff(log(GOOGLE))
adf.test(LOG.GOOGLE.DIFF)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  LOG.GOOGLE.DIFF
## Dickey-Fuller = -2.7132, Lag order = 3, p-value = 0.2865
## alternative hypothesis: stationary

Ngoài ra ta cũng có thể sử dụng một kiểm định khác ngoài ADF test là Phillip-perron (PP test) để kiểm tra tính dừng.

pp.test(GOOGLE.DIFF)

## 
##  Phillips-Perron Unit Root Test
## 
## data:  GOOGLE.DIFF
## Dickey-Fuller Z(alpha) = -53.61, Truncation lag parameter = 3,
## p-value = 0.01
## alternative hypothesis: stationary

Khi kiểm định adf.test() quá chặt để tiêu chuẩn dừng của chuỗi được đáp ứng. Ta có thể sử dụng kiểm định Phillips-Perron. Trong trường hợp này chuỗi thỏa mãn tiêu chuẩn dừng.

2.1. Số dư

Số dư từ phương trình hồi qui ar() có thể thu được thông qua phần tử resid từ mô hình hồi qui vừa khởi tạo.

ar_reg$resid

## Time Series:
## Start = 1 
## End = 60 
## Frequency = 1 
##  [1]          NA          NA          NA  -5.3173676   0.2034711
##  [6]  -7.3559879   1.3485037  19.5382826  -0.3389463  12.7080298
## [11] -12.3518030  11.0575845  15.7948306  12.7202479  -4.9100903
## [16]   5.5718406  -2.1415372   3.8579150 -11.7796720   4.7702872
## [21]  -1.7802753 -55.6748886 -47.2126005  19.7968811 -14.1467656
## [26]  -9.3565211  31.1294524   6.2396480  30.6421896   3.0914063
## [31]  12.5136957   6.3084837   3.6988228   1.8574706  -4.6758812
## [36]  18.3657733   7.8803364 -21.9220692 -13.1731826 -42.1902281
## [41]  29.1869279   8.9884915  23.8531939  10.4428383   9.2800760
## [46]  30.2656575  -2.8932446 -25.7523979   0.6934298  -6.8113620
## [51]  -0.8596133 -34.2914202  -1.4667160   0.3840606 -23.7993729
## [56] -16.6392738  25.2907656 -36.9814898  26.4400337   9.9020792

plot(ar_reg$resid)

Biểu diễn giá trị dự báo và giá trị thực tế

plot(GOOGLE)
points(GOOGLE[30:60]+GOOGLE.DIFF[31:61])

Một trong những giả thuyết quan trọng của mô hình hồi qui chuỗi thời gian đó là phần dư có phân phối chuẩn với kì vọng bằng 0 và không có hiện tượng tự tương quan. Để kiểm định tính phân phối chuẩn của phần dư ta có thể sử dụng kiểm định Jarque-Bera hoặc thông qua đồ thị Q-Q plot.

Kiểm định Jarque-Bera:

Giá trị kiểm định Jarque-Bera được xây dựng dựa trên hệ số skewness và kutorsis theo phương trình: $$JB = \frac{n}{6}(S^2+\frac{(K-3)^2}{4})$$ Trong đó $S = \frac{\hat{\mu_3}}{\hat\mu_2}^{\frac{3}{2}}$ là hệ số skewness và $K = \frac{\hat\mu_4}{\hat\mu_2^2}$ là hệ số kutorsis. Ở đây các chỉ số $\hat\mu_i$ được gọi là moment bậc i của chuỗi được tính theo công thức: $$\hat\mu_i = \frac{1}{n} \sum_{i=1}^{n}(x_i - \bar{x})^i$$ JB sẽ tuân theo qui luật phân phối khi bình phương bởi nó là tổng của bình phương 2 phân phối chuẩn. Miền bác bỏ giả thuyết $H_0$ (chuỗi không có phân phối chuẩn) tại mức level $\alpha$ là $JB \geq \chi_{1-\alpha/2}^2$ hoặc ta có thể sử dụng p-value để so sánh với $\alpha = 0.05$ trong mức ý nghĩa thống kê $1-\alpha = 0.95$ để kết luận bác bỏ $H_0$ nếu p-value nhỏ hơn ngưỡng này. Trong R ta có thể sử dụng kiểm định jarque-bera như sau:

jarque.bera.test(ar_reg$resid[-(1:3)])

## 
##  Jarque Bera Test
## 
## data:  ar_reg$resid[-(1:3)]
## X-squared = 6.4966, df = 2, p-value = 0.03884

hist(ar_reg$resid[-(1:3)], 20)

Ta thấy phân phối của phần dư có đuôi trái dẹp và đuôi phải dày nên nhận định chuỗi không có phân phối chuẩn là khá chính xác. Ngoài ra ta cũng có thể sử dụng biểu đồ Q-Q plot:

qqnorm(ar_reg$resid[-(1:3)])
qqline(ar_reg$resid[-(1:3)])

Đồ thị Q-Q plot sẽ sắp xếp thứ tự các điểm theo ngũ phân vị (Theoretical Quantiles) và sau đó biểu diễn các điểm này trên đồ thị sao cho trục hoành là giá trị phân vị và trục tung là giá trị gốc của chuỗi. Đường chéo chính sẽ là giá trị lý thuyết nếu chuỗi tuân theo phân phối chuẩn. Khi các điểm dữ liệu càng sát với đường chéo chính thì đồ thị càng sát với phân phối chuẩn. Nhìn vào đồ thị Q-Q plot ta có thể thấy phần dư có phân phối đuôi trái quá mỏng (giá trị thực tế nằm dưới giá trị phân phối chuẩn) và đuôi phải dầy. Do đó nó không có dạng phân phối chuẩn.

Kiểm tra yếu tố tự tương quan của phần dư bằng đồ thị ACF thông qua hàm acf():

acf(ar_reg$resid[-(1:3)])

Ta có thể thấy phần dư không có yếu tố tự tương quan do các giá trị tự tương quan không vượt quá ngưỡng tin cậy 95%.

Thông thường các model theo trường phái thống kê sẽ đòi hỏi rất nhiều những tiêu chuẩn về phân phối của phần dư và mối quan hệ giữa các biến. Nổi tiếng trong số những tiêu chuẩn đó là hệ điều kiện Gaussian-Markov để một ước lượng được coi là không chệch tốt nhất (thường được gọi là BLUE - best linear unbiased estimator), bao gồm:

1. Sai số là biến ngẫu nhiên phân phối chuẩn.
2. Kì vọng của sai số bằng 0.
3. Phương sai sai số không đổi.
4. Không có hiện tượng tự tương quan giữa sai số (đối với chuỗi thời gian).
5. Không có hiện tượng đa cộng tuyến giữa các biến dự báo.

Tuy nhiên để thỏa mãn được hệ điều kiện này là rất khó và ngoại trừ những tiêu chuẩn quan trọng như số 2, 5 thì đôi lúc chúng ta chấp nhận một vài tiêu chuẩn không quá quan trọng được phép vi phạm như số 1,3 và 4.

2.2. Dự báo model

Sau khi xây dựng xong model. Để dự báo cho các phiên tiếp theo ta sẽ sử dụng hàm predict(). Tham số cần lưu ý truyền vào bao gồm model hồi qui và số lượng phiên dự báo tiếp theo (n.ahead).

GOOGLE.PRED <- predict(ar_reg, n.ahead = 10)
GOOGLE.PRED

## $pred
## Time Series:
## Start = 61 
## End = 70 
## Frequency = 1 
##  [1] -17.85188520 -12.11731701  15.96896113  -4.30716499 -11.23209622
##  [6]   5.48170080   0.17423671  -7.73212632   0.07579025   0.61788601
## 
## $se
## Time Series:
## Start = 61 
## End = 70 
## Frequency = 1 
##  [1] 19.07721 19.36595 19.69547 21.18716 21.51618 21.69446 21.83206
##  [8] 22.01531 22.05809 22.06036

Kết quả dự báo bao gồm 2 thành phần: pred là giá trị dự báo sai phân bậc nhất của giá chứng khoán và se là sai số chuẩn được dự báo. Để chuyển về dữ chuỗi giá dự báo từ chuỗi sai phân ta sẽ lấy giá trị chứng khoán tại ngày liền cuối cùng cộng với sai phân để thu được ngày dự báo đầu tiên và sử dụng các kết quả của ngày dự báo đầu tiên cộng với sai phân tương ứng để thu được ngày thứ 2,…. Tiếp tục quá trình này đến khi đủ 10 phiên ta sẽ thu được chuỗi giá dự báo.

last_price = as.vector(GOOGLE[length(GOOGLE)])
pred = c(last_price)

for (i in 1:10){
  next_price = pred[length(pred)] + GOOGLE.PRED$pred[i]
  pred = c(pred, next_price)
}
pred

##  [1] 1037.140 1019.288 1007.171 1023.140 1018.833 1007.601 1013.082
##  [8] 1013.256 1005.524 1005.600 1006.218

Đối chiếu với mức giá ở 10 phiên tiếp theo.

actual <- tseries::get.hist.quote(instrument = 'GOOGL', start = '2018-04-02', end = '2018-04-15')

## time series ends   2018-04-13

actual$Close

## 2018-04-02 2018-04-03 2018-04-04 2018-04-05 2018-04-06 2018-04-09 
##    1012.63    1018.68    1029.71    1032.64    1009.95    1020.09 
## 2018-04-10 2018-04-11 2018-04-12 2018-04-13 
##    1036.50    1025.06    1037.29    1036.04

Hiển thị kết quả dự báo và thực tế

plot(as.vector(actual$Close), type = 'l', ylab = 'Value', ylim = c(1000, 1040))
lines(pred[2:11], lty = 'dotted')