Bài học này hướng dẫn bạn cách sử dụng R để thực hiện một số phân tích chuỗi thời gian đơn giản với giả định rằng bạn đã có một số kiến thức cơ bản về R và chuỗi thời gian. Trọng tâm của bài học là hướng dẫn thực hành phân tích chuỗi thời gian trên R chứ không đặt nặng hay đi sâu vào các lý thuyết mô hình chuỗi thời gian. Nội dung tập trung vào các phần: phân rã chuỗi, các phương pháp làm mịn chuỗi, và thực hiện dự báo với mô hình ARIMA. các phương pháp làm mịn chuỗi, và thực hiện dự báo với mô hình ARIMA. Bài viết này được tham khảo và dịch từ Using R for Time Series Analysis.

1 Dữ liệu sử dụng

1.1 Dữ liệu tuổi thọ

Sử dụng hàm scan để đọc dữ liệu dữ liệu dạng .dat vào R. Ở đây, do ba dòng đầu bao gồm các mô tả về dữ liệu, vì vậy, tôi sẽ bỏ qua ba dòng này khi đọc dữ liệu vào R. Sử dụng tùy chọn skip = 3 để bỏ qua ba dòng đầu này.

kings <- scan("http://robjhyndman.com/tsdldata/misc/kings.dat", skip = 3)
kings

##  [1] 60 43 67 50 56 42 50 65 68 43 65 34 47 34 49 41 13 35 53 56 16 43 69 59 48
## [26] 59 86 55 68 51 33 49 67 77 81 67 71 81 68 70 77 56

Thông tin về tuổi thọ của 42 vị vua ở Anh Quốc được đọc và gán vào đối tượng kings trên R. Tiếp đến là khai báo một đối tượng chuỗi thời gian và lưu dữ liệu vào đó. Trên R, chúng ta sử dụng hàm ts như sau:

kingstimeseries <- ts(kings)
kingstimeseries

## Time Series:
## Start = 1 
## End = 42 
## Frequency = 1 
##  [1] 60 43 67 50 56 42 50 65 68 43 65 34 47 34 49 41 13 35 53 56 16 43 69 59 48
## [26] 59 86 55 68 51 33 49 67 77 81 67 71 81 68 70 77 56

1.2 Dữ liệu số ca sinh

Đôi khi tập dữ liệu chuỗi thời gian mà bạn có có thể đã được thu thập theo các tần suất dưới một năm, chẳng hạn như hàng tháng hoặc hàng quý. Chúng ta có thể sử dụng tùy chọn frequency để khai báo tần suất này, chẳng khai báo frequency = 12 cho tần suất thu thập là tháng hoặc frequency = 4 cho các loại dữ liệu quý. Ngoài ra, đối với các chuỗi thời gian dữ liệu được thu thập liên tục thì chúng ta có thể sử dụng tùy chọn start để khai báo thời điểm bắt đầu của dữ liệu phân tích.

Một ví dụ là tập dữ liệu về số ca sinh mỗi tháng ở thành phố New York, từ tháng 1 năm 1946 đến tháng 12 năm 1959 được thu thập lần đầu bởi Newton. Trong trường hợp này, với tần suất thu thập theo tháng, tôi sẽ khai báo ở tùy chọn frequency = 12 và thời điểm bắt đầu phân tích là tháng 1 năm 1946 ở tùy chọn start = c(1946,1)

births <- scan("http://robjhyndman.com/tsdldata/data/nybirths.dat")
birthstimeseries <- ts(births, frequency = 12, start = c(1946, 1))
birthstimeseries

##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 1946 26.663 23.598 26.931 24.740 25.806 24.364 24.477 23.901 23.175 23.227
## 1947 21.439 21.089 23.709 21.669 21.752 20.761 23.479 23.824 23.105 23.110
## 1948 21.937 20.035 23.590 21.672 22.222 22.123 23.950 23.504 22.238 23.142
## 1949 21.548 20.000 22.424 20.615 21.761 22.874 24.104 23.748 23.262 22.907
## 1950 22.604 20.894 24.677 23.673 25.320 23.583 24.671 24.454 24.122 24.252
## 1951 23.287 23.049 25.076 24.037 24.430 24.667 26.451 25.618 25.014 25.110
## 1952 23.798 22.270 24.775 22.646 23.988 24.737 26.276 25.816 25.210 25.199
## 1953 24.364 22.644 25.565 24.062 25.431 24.635 27.009 26.606 26.268 26.462
## 1954 24.657 23.304 26.982 26.199 27.210 26.122 26.706 26.878 26.152 26.379
## 1955 24.990 24.239 26.721 23.475 24.767 26.219 28.361 28.599 27.914 27.784
## 1956 26.217 24.218 27.914 26.975 28.527 27.139 28.982 28.169 28.056 29.136
## 1957 26.589 24.848 27.543 26.896 28.878 27.390 28.065 28.141 29.048 28.484
## 1958 27.132 24.924 28.963 26.589 27.931 28.009 29.229 28.759 28.405 27.945
## 1959 26.076 25.286 27.660 25.951 26.398 25.565 28.865 30.000 29.261 29.012
##         Nov    Dec
## 1946 21.672 21.870
## 1947 21.759 22.073
## 1948 21.059 21.573
## 1949 21.519 22.025
## 1950 22.084 22.991
## 1951 22.964 23.981
## 1952 23.162 24.707
## 1953 25.246 25.180
## 1954 24.712 25.688
## 1955 25.693 26.881
## 1956 26.291 26.987
## 1957 26.634 27.735
## 1958 25.912 26.619
## 1959 26.992 27.897

1.3 Dữ liệu doanh số bán hàng

Tương tự, dữ liệu về doanh số bán hàng hàng tháng của một cửa cửa hàng quà lưu niệm tại một thị trấn nghỉ mát ven biển ở Queensland, Úc, từ tháng 1 năm 1987 đến tháng 12 năm 1993 của Wheelwright and Hyndman, 1998.

souvenir <- scan("http://robjhyndman.com/tsdldata/data/fancy.dat")
souvenirtimeseries <- ts(souvenir, frequency=12, start=c(1987,1))
souvenirtimeseries

##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 1987   1664.81   2397.53   2840.71   3547.29   3752.96   3714.74   4349.61
## 1988   2499.81   5198.24   7225.14   4806.03   5900.88   4951.34   6179.12
## 1989   4717.02   5702.63   9957.58   5304.78   6492.43   6630.80   7349.62
## 1990   5921.10   5814.58  12421.25   6369.77   7609.12   7224.75   8121.22
## 1991   4826.64   6470.23   9638.77   8821.17   8722.37  10209.48  11276.55
## 1992   7615.03   9849.69  14558.40  11587.33   9332.56  13082.09  16732.78
## 1993  10243.24  11266.88  21826.84  17357.33  15997.79  18601.53  26155.15
##            Aug       Sep       Oct       Nov       Dec
## 1987   3566.34   5021.82   6423.48   7600.60  19756.21
## 1988   4752.15   5496.43   5835.10  12600.08  28541.72
## 1989   8176.62   8573.17   9690.50  15151.84  34061.01
## 1990   7979.25   8093.06   8476.70  17914.66  30114.41
## 1991  12552.22  11637.39  13606.89  21822.11  45060.69
## 1992  19888.61  23933.38  25391.35  36024.80  80721.71
## 1993  28586.52  30505.41  30821.33  46634.38 104660.67

2 Vẽ đồ thị chuỗi

Khi bạn đã đọc dữ liệu vào R, bước tiếp theo là tạo các đồ thị mô tả sự biến động hoặc phân phối của chuỗi. Trên R, bạn có thể sử dụng hàm plot.ts() như sau:

plot.ts(kingstimeseries)
plot.ts(birthstimeseries)

Tuổi thọ các vị vua (trái) và số ca sinh (phải)

Chúng ta có thể thấy rằng sự thay đổi trong tuổi thọ của các vị vua (hình bên trái) có thể được mô tả bằng cách sử dụng một mô hình cộng (additive model), vì các dao động ngẫu nhiên trong chuỗi có độ lớn gần như không đổi theo thời gian.

Hoặc ở chuỗi số ca sinh trong tháng (hình bên phải) cho thấy dường như có sự thay đổi theo mùa về số lượng ca sinh mỗi tháng: cao điểm vào mỗi mùa hè và thấp điểm vào mỗi mùa đông. Một lần nữa, có vẻ như chuỗi thời gian này có thể được mô tả bằng cách sử dụng một mô hình cộng, vì mức độ dao động theo mùa gần như không đổi theo thời gian và có vẻ không phụ thuộc vào cấp độ (level) của chuỗi. Các dao động ngẫu nhiên dường như cũng có độ lớn gần như không đổi theo thời gian.

Bên dưới là biểu đồ doanh số hàng tháng của cửa hàng lưu niệm

plot.ts(souvenirtimeseries)
logsouvenirtimeseries <- log(souvenirtimeseries)
plot.ts(logsouvenirtimeseries)

Chuỗi gốc (trái) và chuỗi dạng logarit (phải)

Trong trường hợp này, có vẻ mô hình cộng tính không còn phù hợp để mô tả sự biến động của chuỗi doanh số. Điều này có thể do quy mô của biến động theo mùa và biến động ngẫu nhiên dường như tăng theo cấp độ của chuỗi thời gian. Vì vậy, chúng ta có thể cần phải biến đổi chuỗi thời gian để chuỗi sau khi chuyển đổi là một dạng mô hình cộng tính. Ở đây, tôi sử dụng biến đổi logarit để biến đổi chuỗi ban đầu.

Qua đây có thể thấy độ lớn của các biến động theo mùa và biến động ngẫu nhiên trong chuỗi dạng logarit dường như gần như không đổi theo thời gian và không phụ thuộc trên cấp độ (level) của chuỗi thời gian. Do đó, chuỗi thời gian dạng logarit này có thể được mô tả bằng cách sử dụng một mô hình cộng tính.

3 Phân rã chuỗi

Phân rã một chuỗi thời gian là tách chuỗi thành các cấu phần đặc trưng của chuỗi. Thông thường, một chuỗi thời gian sẽ bao gồm ba thành phần chính là thành phần xu hướng (trend component), thành phần bất thường (irregular component) và thành phần mùa vụ (seasonal component) nếu có.

3.1 Phân rã chuỗi không có yếu tố mùa vụ

Một chuỗi thời gian không có yếu tố mùa vụ (non-seasonal time series) mùa bao gồm một thành phần xu hướng và một thành phần bất thường.

Phân tách chuỗi thời gian liên quan đến việc cố gắng tách chuỗi thời gian thành các các thành phần, nghĩa là ước tính thành phần xu hướng và thành phần bất thường. Để ước tính thành phần xu hướng của chuỗi thời gian không có yếu tố mùa vụ dạng cộng tính, người ta thường sử dụng phương pháp làm mịn (smoothing method), chẳng hạn như tính toán trung bình động đơn giản (simple moving average) của chuỗi.

Hàm SMA() trong gói TTR có thể được sử dụng để làm mượt chuỗi bằng cách sử dụng kỹ thuật trung bình trượt giản đơn. Để sử dụng hàm này, đầu tiên bạn cần cài đặt gói TTR vào máy và gọi thư viện TTR.

library("TTR")

## Warning: package 'TTR' was built under R version 4.2.2

Sau đó, bạn có thể sử dụng hàm SMA() để làm mịn dữ liệu chuỗi thời gian. Đầu tiên là chỉ định thứ tự (span) của đường trung bình động đơn giản, sử dụng tham số “n”. Ví dụ, để tính một đường trung bình động đơn giản bậc 5, chúng ta đặt n = 5 trong hàm SMA().

Trở lại ví dụ về tuổi thọ của 42 vị vua kế tiếp nhau của nước Anh là chuỗi không có yếu tố mùa vụ, được mô tả bằng một mô hình cộng, vì các dao động ngẫu nhiên trong dữ liệu có độ lớn gần như không đổi theo thời gian. Vì vậy, chúng ta có thể ước tính thành phần xu hướng của chuỗi bằng phương pháp làm mịn chuỗi với kỹ thuật trung bình trượt giản đơn. Thiết lập n = 3 để tính toán trung bình 3 kỳ và hiển thị chuỗi đã được làm mịn như sau:

kingstimeseriesSMA3 <- SMA(kingstimeseries, n = 3)
plot.ts(kingstimeseriesSMA3)

Dường như vẫn còn khá nhiều biến động ngẫu nhiên trong chuỗi thời gian được làm mịn với SMA = 3. Do đó, để ước tính thành phần xu hướng chính xác hơn, tôi sẽ thử làm mịn dữ liệu bằng một đường trung bình động đơn giản có bậc cao hơn. Tôi sẽ sai thử vài lần để tìm ra mức độ làm mịn phù hợp. Ví dụ bên dưới là chuỗi bên dưới đã được làm mịn với bậc bằng 8.

kingstimeseriesSMA8 <- SMA(kingstimeseries, n = 8)
plot.ts(kingstimeseriesSMA8)

Dữ liệu được làm mịn với đường trung bình động đơn giản bậc 8 cho một bức tranh rõ ràng hơn về thành phần xu hướng. Chúng ta có thể thấy rằng tuổi thọ của các vị vua Anh dường như đã giảm từ khoảng 55 tuổi xuống còn khoảng 38 tuổi trong các triều đại của 20 vị vua đầu tiên, và sau đó tăng lên khoảng 73 tuổi vào cuối các triều đại của vị vua thứ 40 trong chuỗi thời gian.

3.2 Phân rã chuỗi có yếu tố mùa vụ

Một chuỗi thời gian có yếu tố mùa vụ (seasonal time series) bao gồm thành phần xu hướng, thành phần theo mùa và thành phần bất thường. Vì vậy, phân tách một chuỗi thời có yếu tố mùa vụ sẽ tách chuỗi thời gian thành ba thành phần trên.

Để ước tính thành phần xu hướng và thành phần theo mùa của chuỗi thời gian có yếu tố mùa vụ dạng mô hình cộng, chúng ta có thể sử dụng hàm decompose() trong R. Hàm này ước tính thành phần xu hướng, thành phần mùa vụ và thành phần bất thường của chuỗi.

Kết quả trả về của hàm decompose() sẽ trả về một đối tượng danh sách với kết quả ước tính của các thành phần với tên gọi tương ứng là seasonal, trend, và random.

Ví dụ, như đã thảo luận ở trên, chuỗi thời gian của số ca sinh mỗi tháng ở thành phố New York có yếu tố mùa vụ (nhiều nhất vào mỗi mùa hè và thấp nhất vào mỗi mùa đông) và được mô tả bằng một mô hình cộng (vì các dao động theo mùa và ngẫu nhiên dường như có độ lớn gần như không đổi theo thời gian).

Sử dụng hàm decompose() để phân rã các chuỗi thời gian có yếu tố mùa vụ trên R như sau:

birthstimeseriescomponents <- decompose(birthstimeseries)

Các giá trị ước tính của thành phần mùa vụ, xu hướng và bất thường được lưu trong các biến như birthstimeseriescomponents$seasonal, birthstimeseriescomponents$trend và birthstimeseriescomponents$random.

Ví dụ bên dưới tôi sẽ in ra các giá trị ước tính của thành phần mùa vụ:

birthstimeseriescomponents$seasonal

##             Jan        Feb        Mar        Apr        May        Jun
## 1946 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1947 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1948 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1949 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1950 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1951 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1952 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1953 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1954 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1955 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1956 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1957 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1958 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
## 1959 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
##             Jul        Aug        Sep        Oct        Nov        Dec
## 1946  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1947  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1948  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1949  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1950  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1951  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1952  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1953  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1954  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1955  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1956  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1957  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1958  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197
## 1959  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197

Các yếu tố mùa vụ được ước tính cho các tháng từ tháng 1 đến tháng 12 và giống nhau cho mỗi năm. Tác động mùa vụ lớn nhất là vào tháng 7 (khoảng 1,46) và thấp nhất là vào tháng 2 (khoảng -2,08). Điều này cho thấy dường như có cao điểm về số ca sinh vào tháng 7 và thấp nhất về số ca sinh vào tháng 2 hàng năm.

Chúng ta có thể vẽ đồ thị xu hướng ước tính, theo mùa và các thành phần bất thường của chuỗi thời gian bằng cách sử dụng hàm plot() như sau:

plot(birthstimeseriescomponents)

Biểu đồ trên hiển thị chuỗi thời gian gốc (trên cùng), thành phần xu hướng ước tính (thứ hai từ trên xuống), thành phần mùa vụ ước tính (thứ ba từ trên xuống) và thành phần bất thường ước tính (dưới cùng). Chúng ta thấy rằng thành phần xu hướng ước tính giảm nhẹ từ khoảng 24 vào năm 1947 xuống còn khoảng 22 vào năm 1948. Sau đó là mức tăng ổn định từ đó đến khoảng 27 vào năm 1959.

3.3 Khử yếu tố mùa vụ

Nếu bạn có chuỗi thời gian có yếu tố mùa vụ dạng mô hình cộng thì bạn có thể điều chỉnh chuỗi thời gian theo mùa này bằng cách ước tính thành phần mùa vụ và trừ đi thành phần này ước tính này khỏi chuỗi ban đầu.

Ví dụ: để điều chỉnh theo yếu tố mùa vụ ở chuỗi thời gian của số ca sinh mỗi tháng ở thành phố New York, chúng ta có thể ước tính thành phần mùa vụ này bằng cách sử dụng decompose(), sau đó trừ thành phần ước tính này khỏi chuỗi thời gian ban đầu:

birthstimeseriescomponents <- decompose(birthstimeseries)
birthstimeseriesseasonallyadjusted <- birthstimeseries - birthstimeseriescomponents$seasonal

Sau đó, sử dụng hàm plot() chúng ta có thể hiển thị chuỗi đã điều chỉnh yếu tố mùa vụ như sau:

plot(birthstimeseriesseasonallyadjusted)

Bạn có thể thấy rằng sự biến động theo mùa đã bị loại bỏ khỏi chuỗi thời gian được điều chỉnh theo mùa. Chuỗi thời gian được điều chỉnh theo mùa hiện chỉ chứa thành phần xu hướng và thành phần bất thường.

4 Phương pháp dự báo

Phương pháp làm trơn hàm mũ hay còn gọi là phương pháp làm mịn theo cấp số nhân có thể được sử dụng để đưa ra dự báo ngắn hạn cho dữ liệu chuỗi thời gian.

4.1 Phương pháp làm trơn hàm mũ đơn giản

4.1.1 Giới thiệu

Nếu bạn có một chuỗi thời gian có thể được mô tả bằng mô hình cộng với hằng số và không có tính thời vụ, bạn có thể sử dụng phương pháp làm trơn hàm mũ đơn giản (simple exponential smoothing) để thực hiện các dự báo ngắn hạn.

Phương pháp làm mịn hàm mũ đơn giản cung cấp một cách ước tính mức ở thời điểm hiện tại. Việc làm mịn được điều khiển bởi tham số alpha; để ước tính mức độ tại thời điểm hiện tại. Giá trị của alpha nằm giữa 0 và 1, cho biết tầm quan trọng của các quan sát gần nhất đến giá trị dự báo trong tương lai. Alpha gần bằng 0 có nghĩa là các quan sát gần đây nhất có trọng số nhỏ lên giá trị dự báo.

Ví dụ sử dụng dữ liệu về lượng mưa hàng năm tại Luân Đôn, giai đoạn từ 1813-1912 của Hipel and McLeod, 1994.

rain <- scan("http://robjhyndman.com/tsdldata/hurst/precip1.dat", skip = 1)
rainseries <- ts(rain, start = c(1813))
plot.ts(rainseries)

Từ đồ thị bạn có thể thấy rằng có một mức gần như không đổi (giá trị trung bình vẫn không đổi ở khoảng 25 inch). Các dao động ngẫu nhiên trong chuỗi thời gian dường như có độ lớn gần như không đổi theo thời gian. Vì vậy có lẽ thích hợp để mô tả dữ liệu bằng cách sử dụng một mô hình cộng. Do đó, chúng ta có thể đưa ra dự báo bằng cách sử dụng phương pháp làm mịn hàm mũ đơn giản.

Để đưa ra dự báo bằng cách sử dụng phương pháp làm mịn hàm mũ đơn giản trong R, sử dụng hàm HoltWinters(). Để sử dụng hàm HoltWinters() để làm mịn hàm mũ đơn giản, chúng ta cần đặt tham số beta = FALSE và gamma = FALSE trong hàm. Hàm HoltWinters() trả về một biến danh sách.

4.1.2 Thực hiện

Ví dụ: để sử dụng phương pháp làm mịn hàm mũ đơn giản để đưa ra dự báo về thời gian lượng mưa hàng năm ở Luân Đôn, tôi sử dụng hàm HoltWinters() như sau:

rainseriesforecasts <- HoltWinters(rainseries, beta = FALSE, gamma = FALSE)
rainseriesforecasts

## Holt-Winters exponential smoothing without trend and without seasonal component.
## 
## Call:
## HoltWinters(x = rainseries, beta = FALSE, gamma = FALSE)
## 
## Smoothing parameters:
##  alpha: 0.02412151
##  beta : FALSE
##  gamma: FALSE
## 
## Coefficients:
##       [,1]
## a 24.67819

Kết quả của hàm HoltWinters() cho chúng ta biết rằng giá trị ước tính của tham số alpha là khoảng 0,024. Giá trị này rất gần bằng 0. Điều này cho chúng tôi biết rằng các dự báo dựa trên cả những quan sát gần đây và xa hơn.

Theo mặc định, hàm HoltWinters() chỉ đưa ra dự báo cho cùng khoảng thời gian xem xét. Trong trường hợp này, chuỗi thời gian ban đầu của chúng ta là lượng mưa ở Luân Đôn giai đoạn 1813-1912. Vì vậy các dự báo cũng dành cho giai đoạn 1813-1912.

Trong ví dụ trên, chúng ta đã lưu kết quả của hàm HoltWinters() trong biến danh sách tên là rainseriesforecasts. Các dự báo từ HoltWinters() được lưu trữ là một phần tử có tên fitted trong biến danh sách rainseriesforecasts. Vì vậy chúng ta có thể lấy các giá trị của chúng bằng cách nhập:

head(rainseriesforecasts$fitted, 5)

##          xhat    level
## [1,] 23.56000 23.56000
## [2,] 23.62054 23.62054
## [3,] 23.57808 23.57808
## [4,] 23.76290 23.76290
## [5,] 23.76017 23.76017

Chúng ta có thể vẽ chuỗi thời gian ban đầu được dự báo:

plot(rainseriesforecasts)

Biểu đồ cho thấy chuỗi dự báo (đường màu đỏ) mượt mà hơn nhiều so với chuỗi thời gian gốc (đường màu đen). Để đo độ chính xác của các dự báo, chúng ta có thể tính tổng sai số bình phương (sum of squared errors) cho các sai số dự báo trong mẫu. Nghĩa là các sai số dự báo cho khoảng thời gian được xem xét ở chuỗi thời gian ban đầu của chúng ta. Tổng bình phương các sai số được lưu trữ trong một phần tử được đặt tên SSE của biến danh sách rainseriesforecasts. Chúng ta có thể nhận được giá trị này bằng cách:

rainseriesforecasts$SSE

## [1] 1828.855

Tức là, ở đây tổng bình phương sai số là 1828,855. Việc làm trơn hàm mũ giản đơn thường sử dụng giá trị đầu tiên trong chuỗi thời gian làm giá trị khởi tạo. Ví dụ, trong chuỗi thời gian lượng mưa ở Luân Đôn, giá trị đầu tiên là 23,56 (inch) cho lượng mưa năm 1813. Bạn có thể chỉ định giá trị ban đầu cho hàm HoltWinters() bằng cách sử dụng tham số l.start. Ví dụ, để đưa ra dự báo với giá trị ban đầu là 23,56, chúng ta gõ:

HoltWinters(rainseries, beta = FALSE, gamma = FALSE, l.start = 23.56)

## Holt-Winters exponential smoothing without trend and without seasonal component.
## 
## Call:
## HoltWinters(x = rainseries, beta = FALSE, gamma = FALSE, l.start = 23.56)
## 
## Smoothing parameters:
##  alpha: 0.02412151
##  beta : FALSE
##  gamma: FALSE
## 
## Coefficients:
##       [,1]
## a 24.67819

4.1.3 Dự báo

Như đã giải thích ở trên, theo mặc định, hàm HoltWinters() chỉ đưa ra dự báo cho khoảng thời gian được bao quát của dữ liệu gốc. Đó là giai đoạn 1813-1912 cho chuỗi thời gian lượng mưa. Chúng ta có thể đưa ra dự báo cho các mốc thời gian tiếp theo bằng cách sử dụng hàm forecast.HoltWinters() trong gói forecast trên R.

library("forecast")

## Warning: package 'forecast' was built under R version 4.2.2

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

Đối số đầu tiên (đầu vào) của hàm forecast.HoltWinters() là kết quả dự báo từ hàm HoltWinters().

Ví dụ: trong trường hợp chuỗi thời gian lượng mưa, chúng ta đã lưu trữ mô hình dự báo bằng hàm HoltWinters() ở biến rainseriesforecasts. Bạn chỉ định số lượng điểm thời gian tiếp theo mà bạn muốn đưa ra dự báo bằng cách sử dụng tham số h trong hàm forecast.HoltWinters(). Bên dưới là phần dự báo lượng mưa cho các năm 1814-1820 (8 năm tiếp theo) bằng cách sử dụng forecast.HoltWinters():

rainseriesforecasts2 <- forecast:::forecast.HoltWinters(rainseriesforecasts, h = 8) 
rainseriesforecasts2

##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1913       24.67819 19.17493 30.18145 16.26169 33.09470
## 1914       24.67819 19.17333 30.18305 16.25924 33.09715
## 1915       24.67819 19.17173 30.18465 16.25679 33.09960
## 1916       24.67819 19.17013 30.18625 16.25434 33.10204
## 1917       24.67819 19.16853 30.18785 16.25190 33.10449
## 1918       24.67819 19.16694 30.18945 16.24945 33.10694
## 1919       24.67819 19.16534 30.19105 16.24701 33.10938
## 1920       24.67819 19.16374 30.19265 16.24456 33.11182

Hàm forecast.HoltWinters() cung cấp cho bạn dự báo cho một năm với khoảng dự báo 80% và 95%. Ví dụ, lượng mưa dự báo cho năm 1920 là khoảng 24,68 inch, với khoảng dự báo 95% là (16.24, 33.11).

Để vẽ các dự báo được đưa ra bởi forecast.HoltWinters(), chúng ta có thể sử dụng hàm plot.forecast():

forecast:::plot.forecast(rainseriesforecasts2)

Ở đây, các dự báo cho giai đoạn 1913-1920 được vẽ dưới dạng một đường màu xanh lam, khoảng dự báo 80% dưới dạng vùng được tô bóng màu cam và khoảng dự báo 95% là vùng được tô bóng màu vàng. Các sai số dự báo được tính bằng giá trị quan sát được trừ đi giá trị dự báo cho mỗi thời điểm. Chúng tôi chỉ có thể tính toán các sai số dự báo trong khoảng thời gian được đề cập theo chuỗi thời gian ban đầu.

4.1.4 Kiểm định

Như được đề cập ở trên, một phép đo độ chính xác của mô hình dự báo là tổng bình phương sai số (SSE) cho các sai số dự báo trong mẫu (in-sample forecast errors).

Các sai số dự báo trong mẫu được lưu ở một phần tử có tên là residuals của biến danh sách được trả về bởi forecast.HoltWinters(). Nếu mô hình dự báo không thể được cải thiện, không nên có mối tương quan giữa các sai số dự báo cho các dự báo liên tiếp.

Nói cách khác, nếu có mối tương quan giữa các sai số dự báo đối với các dự báo liên tiếp, có khả năng là các dự báo làm mịn hàm mũ đơn giản có thể được cải thiện bởi một kỹ thuật dự báo khác.

4.1.4.1 Sự tự tương quan

Để tìm hiểu sự tương quan này, chúng ta có thể sử dụng biểu đồ tương quan (correlogram) của các sai số trong mẫu cho các độ trễ 1-20. Chúng ta có thể tính toán biểu đồ tương quan của các sai số dự báo bằng cách sử dụng hàm acf() trên R. Để chỉ định độ trễ tối đa mà chúng tôi muốn xem xét, chúng tôi sử dụng tùy chọn lag.max.

Ví dụ: để tính biểu đồ tương quan của các sai số dự báo trong mẫu cho chuỗi lượng mưa ở Luân Đôn với độ trễ 1-20:

acf(rainseriesforecasts2$residuals, lag.max = 20, na.action = na.pass)

Bạn có thể thấy từ biểu đồ tương quan mẫu rằng tự tương quan ở độ trễ 3 chỉ chạm vào các giới hạn ý nghĩa thống kê. Để kiểm tra xem có bằng chứng quan trọng nào cho các mối tương quan khác không hay không ở độ trễ 1-20, chúng ta có thể thực hiện kiểm định Ljung-Box. Kiểm định Ljung-Box được thực hiện trong R qua hàm Box.test(). Độ trễ tối đa mà chúng ta muốn xem xét được thiết lập ở tùy chọn lag.

Ví dụ bên dưới sẽ kiểm tra xem có tự tương quan khác không tại các độ trễ 1-20 cho với các sai số dự báo trong mẫu của chuỗi liệu lượng mưa ở Luân Đôn.

Box.test(rainseriesforecasts2$residuals, lag = 20, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  rainseriesforecasts2$residuals
## X-squared = 17.401, df = 20, p-value = 0.6268

Kết quả giá trị thống kê Ljung-Box là 17,4 và giá trị p là 0,6. Điều đó cho thấy có rất ít bằng chứng về tự tương quan khác không trong các sai số dự báo trong mẫu ở độ trễ 1-20.

Để chắc chắn rằng mô hình dự báo không thể cải thiện được nữa, bạn cũng nên kiểm tra liệu các sai số dự báo có phân phối chuẩn với giá trị trung bình bằng 0 và phương sai không đổi hay không.

4.1.4.2 Tính chất phân phối chuẩn

Để kiểm tra xem các sai số dự báo có phương sai không đổi hay không, chúng ta có đánh giá qua biểu đồ thời gian của các sai số dự báo trong mẫu.

plot.ts(rainseriesforecasts2$residuals)

Biểu đồ cho thấy rằng các sai số dự báo trong mẫu dường như có phương sai gần như không đổi theo thời gian, mặc dù độ lớn của các biến động khi bắt đầu chuỗi thời gian (1820-1830) có thể hơi ít hơn vào những thời điểm sau đó (chẳng hạn, 1840-1850).

Để kiểm tra xem các sai số dự báo có được phân phối chuẩn với giá trị trung bình bằng 0 hay không, chúng ta có thể vẽ biểu đồ của các sai số dự báo, với một đường cong phân phối chuẩn được vẽ chồng lên.

Để làm điều này, chúng ta có thể định nghĩa một hàm R plotForecastErrors() như sau:

plotForecastErrors <- function(forecasterrors)
     {
         # make a histogram of the forecast errors: 
         mybinsize <- IQR(forecasterrors, na.rm = TRUE)/4
         mysd   <- sd(forecasterrors, na.rm = TRUE)
         mymin  <- min(forecasterrors, na.rm = TRUE) - mysd*5      
         mymax  <- max(forecasterrors, na.rm = TRUE) + mysd*3     
         # generate normally distributed data with mean 0 and standard deviation mysd
         mynorm <- rnorm(10000, mean=0, sd = mysd)
         mymin2 <- min(mynorm)
         mymax2 <- max(mynorm)
         if (mymin2 < mymin) { mymin <- mymin2 }
         if (mymax2 > mymax) { mymax <- mymax2 }
         # make a red histogram of the forecast errors, with the normally distributed data overlaid: 
         mybins <- seq(mymin, mymax, mybinsize)
         hist(forecasterrors, col = "red", freq = FALSE, breaks = mybins) 
         # freq=FALSE ensures the area under the histogram = 1
         # generate normally distributed data with mean 0 and standard deviation mysd
         myhist <- hist(mynorm, plot = FALSE, breaks = mybins) 
         # plot the normal curve as a blue line on top of the histogram of forecast errors:
         points(myhist$mids, myhist$density, type = "l", col = "blue", lwd = 2) 
     }

Bạn sẽ sử dụng plotForecastErrors() để vẽ biểu đồ histogram với đường cong phân phối chuẩn cho các sai số dự báo của chuỗi lượng mưa như sau:

plotForecastErrors(rainseriesforecasts2$residuals)

Biểu đồ cho thấy rằng phân phối của các sai số dự báo gần như tập trung vào 0 và ít nhiều có phân phối chuẩn, mặc dù nó có vẻ hơi lệch về bên phải so với một đường cong phân phối chuẩn. Tuy nhiên, độ lệch bên phải là tương đối nhỏ, và do đó có thể xem các sai số dự báo có phân phối xấp xỉ phân phối bình chuẩn với giá trị trung bình bằng 0.

Kiểm định Ljung-Box cho thấy có rất ít bằng chứng về tự tương quan khác 0 trong các sai số mẫu dự báo và phân phối của sai số dự báo xấp xỉ phân phối chuẩn với giá trị trung bình bằng 0. Điều này cho thấy rằng phương pháp làm mịn hàm mũ đơn giản cung cấp một mô hình dự báo đầy đủ cho lượng mưa ở Luân Đôn. Vì vậy, mô hình dự báo có thể là tối ưu hay không thể cải thiện thêm. Hơn nữa, các giả định về khoảng dự báo 80% và 95% dựa trên giả định phân phối chuẩn đều hợp lệ.

4.2 Phương pháp làm mịn hàm mũ Holt

4.2.1 Giới thiệu

Nếu bạn có một chuỗi thời gian có thể được mô tả bằng mô hình cộng với thành phần xu hướng và không có tính thời vụ, bạn có thể sử dụng phương pháp làm mịn theo cấp số nhân của Holt để tạo ra các dự báo ngắn hạn.

Phương pháp làm mịn theo cấp số nhân của Holt ước tính level và độ dốc tại thời điểm hiện tại. Phương pháp làm mịn này được kiểm soát bởi hai tham số, alpha, để ước tính level tại thời điểm hiện tại, và beta để ước tính độ dốc b của thành phần xu hướng tại thời điểm hiện tại.

Giống như phương pháp làm mịn hàm mũ đơn giản, các tham số alpha và beta có các giá trị từ 0 đến 1 với ý nghĩa các giá trị gần 0 thì quan sát gần nhất càng ít quan trọng (trọng số nhỏ) đối với giá trị được dự báo.

4.2.2 Thực hiện

Một ví dụ về chuỗi thời gian có thể được mô tả bằng mô hình cộng với yếu tố xu hướng (nhưng không có yếu tố mùa vụ) là dữ liệu hàng năm của đường kính viền áo phụ nữ tại hem, từ 1866 đến 1911. Dữ liệu có sẵn trong tệp Hipel and McLeod, 1994.

skirts <- scan("http://robjhyndman.com/tsdldata/roberts/skirts.dat", skip = 5)
skirtsseries <- ts(skirts, start = c(1866))
plot.ts(skirtsseries)

Chúng ta có thể thấy từ biểu đồ rằng đã có sự gia tăng đường kính viền áo từ khoảng 600 trong năm 1866 đến khoảng 1050 vào năm 1880, và sau đó đường kính này giảm xuống còn khoảng 520 trong năm 1911.

Để đưa ra dự báo, chúng ta có thể điều chỉnh mô hình dự báo bằng cách sử dụng hàm HoltWinters() trong R. Để sử dụng HoltWinters() cho việc làm trơn hàm mũ của Holt, chúng ta cần đặt tham số gamma=FALSE.

Ví dụ sau sẽ sử dụng phương pháp làm mịn theo cấp số nhân của Holt để ước lượng mô hình dự báo đường kính đường viền.

skirtsseriesforecasts <- HoltWinters(skirtsseries, gamma = FALSE)
skirtsseriesforecasts

## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = skirtsseries, gamma = FALSE)
## 
## Smoothing parameters:
##  alpha: 0.8383481
##  beta : 1
##  gamma: FALSE
## 
## Coefficients:
##         [,1]
## a 529.308585
## b   5.690464

skirtsseriesforecasts$SSE

## [1] 16954.18

Giá trị ước tính của alpha là 0,84 và của beta là 1,00. Cả hai đều cao. Điều đó cho chúng tôi biết rằng cả ước tính giá trị hiện tại cho level và độ dốc b của thành phần xu hướng, chủ yếu dựa trên các quan sát rất gần đây. Điều này có ý nghĩa trực quan tốt, vì level và độ dốc của chuỗi thời gian đều thay đổi khá nhiều theo thời gian. Các giá trị của tổng các sai số bình phương đối với các sai số dự báo trong mẫu là 16954.

Chúng ta có thể vẽ chuỗi thời gian gốc dưới dạng một đường màu đen, với các giá trị dự báo là một đường màu đỏ, bằng cách:

plot(skirtsseriesforecasts)

Từ đồ thị, chúng ta có thể thấy rằng các dự báo trong mẫu phù hợp khá tốt với các giá trị quan sát được, mặc dù chúng có xu hướng tụt lại phía sau các giá trị quan sát được một chút. Nếu muốn, bạn có thể chỉ định các giá trị ban đầu của level và độ dốc b của thành phần xu hướng bằng cách sử dụng tùy chọn l.start và b.start ở hàm HoltWinters(). Người ta thường đặt giá trị ban đầu của level bằng giá trị đầu tiên trong chuỗi thời gian và giá trị ban đầu của độ dốc thành giá trị thứ hai trừ đi giá trị đầu tiên.

Ví dụ, để ước lượng một mô hình dự báo với dữ liệu đường kính viền bằng cách sử dụng phương pháp làm mịn hàm mũ của Holt, với các giá trị ban đầu của 608 cho level và 9 cho độ dốc b của thành phần xu hướng, chúng ta gõ:

HoltWinters(skirtsseries, gamma = FALSE, l.start = 608, b.start = 9)

## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = skirtsseries, gamma = FALSE, l.start = 608, b.start = 9)
## 
## Smoothing parameters:
##  alpha: 0.8346775
##  beta : 1
##  gamma: FALSE
## 
## Coefficients:
##         [,1]
## a 529.278637
## b   5.670129

4.2.3 Dự báo

Đối với việc làm mịn theo cấp số nhân đơn giản, chúng ta có thể đưa ra dự báo cho thời gian trong tương lai ngoài phạm vi của chuỗi thời gian ban đầu bằng cách sử dụng hàm forecast.HoltWinters().

Ví dụ: chúng ta có thể dự báo đường kính viền váy cho giai đoạn từ năm 1912 đến 1930 dựa trên dữ liệu chuỗi ban đầu chúng ta có từ năm 1866 đến năm 1911 như sau:

skirtsseriesforecasts2 <- forecast:::forecast.HoltWinters(skirtsseriesforecasts, h = 19)
forecast:::plot.forecast(skirtsseriesforecasts2)

Đường dự báo được hiển thị dưới dạng đường màu xanh lam, với khoảng dự báo 80% là màu cam khu vực bóng mờ và khoảng dự báo 95% dưới dạng khu vực bóng mờ màu vàng.

4.2.4 Kiểm định

Tương tự với phương pháp làm mịn hàm mũ đơn giản, chúng ta có thể kiểm tra xem mô hình dự báo có thể được cải thiện bằng cách kiểm tra xem các sai số dự báo trong mẫu có sự tự tương quan khác 0 hay không ở các độ trễ 1-20.

4.2.4.1 Sự tự tương quan

Ví dụ: đối với dữ liệu đường viền váy, chúng ta có thể tạo biểu đồ tương quan và thực hiện kiểm định Ljung-Box như sau:

acf(skirtsseriesforecasts2$residuals, lag.max = 20, na.action = na.pass)

Box.test(skirtsseriesforecasts2$residuals, lag = 20, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  skirtsseriesforecasts2$residuals
## X-squared = 19.731, df = 20, p-value = 0.4749

Ở đây biểu đồ tương quan cho thấy tự tương quan mẫu đối với các sai số dự báo trong mẫu ở độ trễ 5 vượt quá các giới hạn ý nghĩa thống kê. Tuy nhiên, chúng tôi kì vọng một trong 20 trường hợp tự tương quan đối với hai mươi độ trễ đầu tiên vượt quá giới hạn ý nghĩa 95% một cách tình cờ. Thật vậy, kết quả kiểm định Ljung-Box với giá trị p là 0,47, cho thấy có rất ít bằng chứng về giá trị tự tương quan khác 0 trong các sai số dự báo trong mẫu ở độ trễ 1-20.

4.2.4.2 Tính chất phân phối chuẩn

Tương tự với phương pháp làm mịn hàm mũ đơn giản, chúng ta cũng nên kiểm tra xem phương sai sai số dự báo có không đổi theo thời gian và có phân phối chuẩn với giá trị trung bình bằng 0 hay không? Chúng ta có thể làm điều này bằng cách đánh giá biểu đồ sai số dự báo và biểu đồ phân phối sai số dự báo bằng hàm plotForecastErrors đã viết ở trên.

plot.ts(skirtsseriesforecasts2$residuals)

plotForecastErrors(skirtsseriesforecasts2$residuals)

Biểu đồ thời gian của các sai số dự báo cho thấy các sai số dự báo có phương sai gần như không đổi theo thời gian. Biểu đồ phân phối của sai số dự báo cho thấy có thể sai số dự báo có phân phối chuẩn với giá trị trung bình bằng 0 và phương sai không đổi.

Do đó, kiểm định Ljung-Box cho thấy có rất ít bằng chứng về hiện tượng tự tương quan trong các sai số dự báo, trong khi biểu đồ thời gian và biểu đồ phân phối của các sai số dự báo cho thấy rằng có thể các sai số dự báo có phân phối xấp xỉ phân phối chuẩn với giá trị trung bình bằng 0 và phương sai không đổi. Ngoài ra, các giả định rằng khoảng dự đoán 80% và 95% được đều hợp lệ. Do đó, chúng ta có thể kết luận rằng phương pháp làm mịn theo cấp số nhân của Holt cung cấp một mô hình dự đoán đầy đủ cho đường kính váy.

4.3 Phương pháp làm trơn hàm mũ Holt-Winters

4.3.1 Giới thiệu

Nếu bạn có một chuỗi thời gian có thể được mô tả bằng mô hình cộng với thành phần xu hướng và có yếu tố mùa vụ, bạn có thể sử dụng phương pháp làm mịn hàm mũ Holt-Winters để đưa ra các dự báo ngắn hạn.

Phương pháp làm mịn hàm mũ Holt-Winters ước tính level, độ dốc và thành phần mùa tại thời điểm hiện tại.

Việc làm mịn được kiểm soát bởi ba tham số: alpha, beta và gamma, để ước tính level, độ dốc b của thành phần xu hướng và thành phần mùa tương ứng tại thời điểm hiện tại. Những thông số alpha, beta và gamma đều có các giá trị từ 0 đến 1 và các giá trị càng gần 0 cho biết các quan sát gần nhất càng ít quan trọng (trọng số nhỏ) đối với các giá trị dự báo trong tương lai.

4.3.2 Thực hiện

Một ví dụ về chuỗi thời gian có thể được mô tả bằng mô hình cộng có xu hướng và tính thời vụ là chuỗi thời gian doanh thu hàng tháng của cửa hàng lưu niệm tại một thị trấn nghỉ mát ven biển ở Queensland, Úc

Để đưa ra dự báo, chúng ta có thể điều chỉnh mô hình dự đoán bằng cách sử dụng hàm HoltWinters().

Ví dụ bên dưới sẽ thực hiện ước lượng mô hình dự báo cho chuỗi doanh thu bán hàng hàng tháng trong cửa hàng lưu niệm:

logsouvenirtimeseries <- log(souvenirtimeseries)
souvenirtimeseriesforecasts <- HoltWinters(logsouvenirtimeseries)
souvenirtimeseriesforecasts

## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = logsouvenirtimeseries)
## 
## Smoothing parameters:
##  alpha: 0.413418
##  beta : 0
##  gamma: 0.9561275
## 
## Coefficients:
##            [,1]
## a   10.37661961
## b    0.02996319
## s1  -0.80952063
## s2  -0.60576477
## s3   0.01103238
## s4  -0.24160551
## s5  -0.35933517
## s6  -0.18076683
## s7   0.07788605
## s8   0.10147055
## s9   0.09649353
## s10  0.05197826
## s11  0.41793637
## s12  1.18088423

souvenirtimeseriesforecasts$SSE

## [1] 2.011491

Các giá trị ước tính của alpha, beta và gamma lần lượt là 0,41, 0,00 và 0,96. Các giá trị của alpha (0,41) tương đối thấp, cho thấy ước tính của mức tại thời điểm hiện tại điểm dựa trên cả những quan sát gần đây và một số quan sát trước đó. Giá trị của beta là 0,00, cho biết rằng ước tính độ dốc b của thành phần xu hướng không được cập nhật theo chuỗi thời gian và thay vào đó được đặt bằng giá trị ban đầu của nó. Điều này có ý nghĩa trực quan tốt, vì level thay đổi khá nhiều trong chuỗi thời gian, nhưng độ dốc b của thành phần xu hướng vẫn gần như giống nhau. Ngược lại, giá trị của gamma (0,96) cao cho thấy ước tính của thành phần mùa ở thời điểm hiện tại chỉ dựa trên những quan sát gần đây nhất.

Tương tự như phương pháp làm trơn hàm mũ đơn giản hay làm trơn hàm mũ Holt, chúng ta có thể vẽ chuỗi thời gian ban đầu dưới dạng một đường màu đen, với các giá trị được dự báo là một đường màu đỏ ở trên cùng:

plot(souvenirtimeseriesforecasts)

Từ biểu đồ, chúng ta thấy rằng phương pháp mũ Holt-Winters rất thành công trong việc dự đoán cao điểm theo mùa, xảy ra vào khoảng tháng 11 hàng năm.

4.3.3 Dự báo

Để đưa ra dự báo cho các thời điểm trong tương lai không có trong chuỗi thời gian ban đầu, chúng ta sử dụng hàm forecast.HoltWinters().

Ví dụ sau sẽ thực hiện dự báo doanh số bán hàng lưu niệm cho tháng 1 năm 1994 đến tháng 12 năm 1998 (48 tháng nữa) và vẽ sơ đồ dự báo dựa trên dữ liệu được thu thập ban đầu (từ tháng 1 năm 1987 đến tháng 12 năm 1993).

souvenirtimeseriesforecasts2 <- forecast:::forecast.HoltWinters(souvenirtimeseriesforecasts, h=48)
forecast:::plot.forecast(souvenirtimeseriesforecasts2)

Các dự báo được hiển thị dưới dạng đường màu xanh lam và các khu vực được tô bóng màu cam và màu vàng hiển thị 80% và khoảng dự đoán 95%, tương ứng.

4.3.4 Kiểm định

4.3.4.1 Sự tự tương quan

Chúng ta có thể đánh giá liệu mô hình dự đoán có thể được cải thiện hay không bằng cách kiểm tra xem liệu sai số dự báo trong mẫu có sự tự tương quan khác không ở độ trễ 1-20, bằng cách vẽ biểu đồ tương quan và thực hiện kiểm định Ljung-Box:

acf(souvenirtimeseriesforecasts2$residuals, lag.max = 20, na.action = na.pass)

Box.test(souvenirtimeseriesforecasts2$residuals, lag = 20, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  souvenirtimeseriesforecasts2$residuals
## X-squared = 17.53, df = 20, p-value = 0.6183

Biểu đồ tương quan cho thấy rằng hệ số tự tương quan đối với các sai số dự báo trong mẫu không vượt quá giới hạn ý nghĩa cho độ trễ 1-20. Hơn nữa, giá trị p của kiểm định Ljung-Box là 0,6, cho thấy có rất ít bằng chứng về tự tương quan khác không ở độ trễ 1-20.

4.3.4.2 Tính chất phân phối chuẩn

Chúng ta có thể kiểm tra xem các sai số dự báo có phương sai không đổi theo thời gian và có phân phối chuẩn với giá trị trung bình bằng 0 hay không? Bằng cách tạo biểu đồ thời gian của các sai số dự báo và biểu đồ phân phối chuẩn với hàm plotForecastErrors đã viết ở trên.

plot.ts(souvenirtimeseriesforecasts2$residuals)

plotForecastErrors(souvenirtimeseriesforecasts2$residuals)

Từ biểu đồ thời gian, có vẻ hợp lý khi các sai số dự báo có phương sai không đổi theo thời gian.

Từ biểu đồ phân phối cho thấy các sai số dự báo có phân phối xấp xỉ phân phối chuẩn với giá trị trung bình bằng 0.

Do đó, có rất ít bằng chứng về sự tự tương quan ở độ trễ 1-20 đối với các sai số dự báo, và các sai số dự báo có phân phối xấp xỉ phân phối chuẩn với giá trị trung bình bằng 0 và phương sai không đổi theo thời gian.

Điều này gợi ý rằng phương pháp làm mịn hàm mũ Holt-Winters cung cấp một mô hình dự báo đầy đủ về doanh số bán hàng tại cửa hàng lưu niệm, mà có lẽ không thể cải thiện thêm được. Hơn nữa, các giả định dựa trên các khoảng dự đoán 80% hay 95% đều hợp lệ.

5 Mô hình ARIMA

5.1 Giới thiệu mô hình ARIMA

Các phương pháp làm mịn theo cấp số nhân rất hữu ích để đưa ra dự báo và không đưa ra giả định nào về mối tương quan giữa các giá trị liên tiếp của chuỗi thời gian. Tuy nhiên, nếu bạn muốn thực hiện các dự báo được thực hiện bằng phương pháp làm mịn theo cấp số nhân, việc dự báo sẽ yêu cầu các sai số dự báo không có sự tương quan và có phân phối chuẩn với trung bình bằng 0 và phương sai không đổi.

Trong khi các phương pháp làm mịn hàm mũ không đưa ra bất kỳ giả định nào về mối tương quan giữa các giá trị của chuỗi thời gian, trong một số trường hợp, bạn có thể tạo mô hình dự báo tốt hơn bằng cách xét tới các mối tương quan trong dữ liệu. Các mô hình tự hồi quy trung bình trượt tích hợp (ARIMA - Autoregressive Integrated Moving Average) là một dạng mô hình thống kê xét tới thành phần bất thường của chuỗi (irregular component), cho phép tồn tại sự tự tương quan trong thành phần bất thường.

5.2 Lấy sai phân một chuỗi thời gian

Các mô hình ARIMA được xác định cho các chuỗi thời gian dừng. Do đó, nếu bạn bắt đầu với một chuỗi không dừng (non-stationary) thì trước tiên bạn sẽ cần lấy sai phân chuỗi cho đến khi bạn có được chuỗi dừng. Nếu bạn phải lấy sai phân chuỗi d lần để có được chuỗi dừng, thì mô hình ARIMA sẽ có dạng tổng quát là ARIMA(p,d,q). Trong đó, d được gọi là bậc tích hợp. Trong thực tế, thông thường chúng ta chỉ lấy sai phân tối đa là đến bậc 2.

Trên R, bạn có thể sử dụng hàm diff() để lấy sai phân của một chuỗi. Ví dụ: chuỗi thời gian của đường kính viền váy phụ nữ, từ 1866 đến 1911 là không dừng, vì level thay đổi rất nhiều theo thời gian:

skirtsseriesdiff1 <- diff(skirtsseries, differences = 1)
plot.ts(skirtsseriesdiff1)

Chuỗi sai phân bậc 1 (ở trên) dường như không dừng. Chúng ta tiếp tục kiểm tra tính dừng của chuỗi sai phân bậc 2:

skirtsseriesdiff2 <- diff(skirtsseries, differences = 2)
plot.ts(skirtsseriesdiff2)

5.3 Kiểm định tính dừng

Các kiểm định chính thức về tính dừng được gọi là kiểm tra nghiệm đơn vị (unit root tests) có trong gói fUnitRoots. Chuỗi sai phân bậc 2 ở trên là dừng theo giá trị trung bình và phương sai, vì level của chuỗi gần như không đổi theo thời gian và phương sai của chuỗi xuất hiện gần như không đổi theo thời gian. Trong trường hợp này, bạn có thể sử dụng mô hình ARIMA(p,2,q) cho chuỗi thời gian trên. Bước tiếp theo là xác định các giá trị p và q cho mô hình ARIMA.

Một ví dụ khác là chuỗi thời gian về tuổi qua đời của các vị vua kế vị của nước Anh:

kingtimeseriesdiff1 <- diff(kingstimeseries, differences = 1)
plot.ts(kingtimeseriesdiff1)

Chuỗi sai phân bậc 1 có vẻ là dừng theo giá trị trung bình và phương sai. Do đó, mô hình ARIMA(p,1,q) có lẽ phù hợp với chuỗi thời gian về tuổi qua đời của các vị vua nước Anh.

Bằng cách lấy sai phân bậc 1 của chuỗi, chúng ta đã loại bỏ thành phần xu hướng của chuỗi và chỉ còn lại thành phần bất thường.

Bây giờ chúng ta có thể kiểm tra xem có mối tương quan nào giữa các số hạng liên tiếp của thành phần bất thường này hay không. Nếu có, điều này có thể giúp chúng ta tạo ra một mô hình dự đoán về tuổi qua đời của các vị vua.

5.4 Lựa chọn mô hình ARIMA

Nếu chuỗi thời gian của bạn là dừng (bậc gốc hoặc sai phân), bước tiếp theo là chọn mô hình ARIMA thích hợp. Điều này nghĩa là tìm giá trị của các giá trị phù hợp nhất của p và q cho mô hình ARIMA(p,d,q). Thông thường, để làm điều này bạn xem xét biểu đồ tương quan (correlogram) và biểu đồ tương quan từng phần (partial correlogram) của các chuỗi dừng.

Để vẽ biểu đồ tương quan và biểu đồ tương quan từng phần, chúng ta lần lượt sử dụng các hàm acf() và hàm pacf() trong R.

Trở lại ví dụ về tuổi thọ của các vị vua ở Anh, để vẽ biểu đồ tương quan cho 20 độ trễ đầu tiên của chuỗi sai phân bậc 1, tôi thực hiện như sau:

acf(kingtimeseriesdiff1, lag.max = 20, na.action = na.pass)

acf(kingtimeseriesdiff1, lag.max = 20, plot = FALSE)

## 
## Autocorrelations of series 'kingtimeseriesdiff1', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000 -0.360 -0.162 -0.050  0.227 -0.042 -0.181  0.095  0.064 -0.116 -0.071 
##     11     12     13     14     15     16     17     18     19     20 
##  0.206 -0.017 -0.212  0.130  0.114 -0.009 -0.192  0.072  0.113 -0.093

Từ biểu đồ tương quan, chúng tôi thấy rằng tự tương quan ở độ trễ 1 (-0,360) vượt quá giới hạn ý nghĩa thống kê nhưng tất cả các tự tương quan khác giữa độ trễ 1-20 không vượt quá giới hạn ý nghĩa thống kê.

Để vẽ biểu đồ tương quan từng phần với 20 độ trễ đầu tiên cho chuỗi sai phân bậc 1 tính các giá trị của tự tương quan từng phần, tôi sử dụng hàm pacf như sau:

pacf(kingtimeseriesdiff1, lag.max = 20, na.action = na.pass)

pacf(kingtimeseriesdiff1, lag.max = 20, plot=FALSE)

## 
## Partial autocorrelations of series 'kingtimeseriesdiff1', by lag
## 
##      1      2      3      4      5      6      7      8      9     10     11 
## -0.360 -0.335 -0.321  0.005  0.025 -0.144 -0.022 -0.007 -0.143 -0.167  0.065 
##     12     13     14     15     16     17     18     19     20 
##  0.034 -0.161  0.036  0.066  0.081 -0.005 -0.027 -0.006 -0.037

Biểu đồ tương quan từng phần cho thấy tự tương quan từng phần ở độ trễ 1, 2 và 3 vượt quá giới hạn ý nghĩa thống kê, là âm, và đang giảm dần về độ lớn theo sự gia tăng ở các độ trễ (độ trễ 1: -0,360, độ trễ 2: -0,335, độ trễ 3: -0,321). Tự tương quan từng phần giảm dần về 0 sau độ trễ 3.

Vì biểu đồ tương quan bằng 0 sau độ trễ 1 và biểu đồ tương quan từng phần giảm dần về 0 sau độ trễ 3. Điều này có nghĩa là các mô hình ARIMA sau đây có thể được sử dụng cho chuỗi sai phân bậc 1 đang xét:

Mô hình ARIMA(3,1,0) với thành phần tự hồi quy bậc 3 (p = 3), vì một biểu đồ tự tương quan từng phần bằng 0 sau độ trễ 3 và biểu đồ tự tương quan kết thúc bằng 0.
Mô hình ARIMA(0,1,1) với thành phần trung bình trượt bậc 1 (q = 1), vì biểu đồ tự tương quan bằng 0 sau độ trễ đầu tiên và biểu đồ tự tương quan từng phần tắt dần về 0.
Mô hình ARIMA(p,1,q), tức là một mô hình hỗn hợp với p và q lớn hơn 0, do biểu đồ tự tương quan và biểu đồ tương quan từng phần giảm dần về 0.

Chúng ta sử dụng nguyên tắc tối giản (parsimony) để quyết định chọn mô hình. Cụ thể, trong số các mô hình phù hợp thì mô hình đơn giản nhất (ít tham số nhất) sẽ là mô hình được chọn.

Ví dụ, mô hình ARMA(3,0) có 3 tham số, mô hình ARMA(0,1) có 1 tham số và mô hình ARMA(p,q) có ít nhất 2 tham số thì mô hình được chọn sẽ là là mô hình trung bình động bậc 1 hay MA(1). Mô hình này có thể được viết lại như sau:

\[X_t - \mu = Z_t - (\theta*Z_{t-1})\]

trong đó: - $X_t$ là một chuỗi dừng - $\mu$ là giá trị trung bình của chuỗi $X_t$, - $Z_t$ là một nhiễu trắng (white noise) với trung bình bằng 0 và phương sai không đổi, - $\theta$ là một tham số cần được ước lượng.

Một mô hình MA (trung bình trượt) được sử dụng để lập mô hình chuỗi thời gian cho thấy sự phụ thuộc ngắn hạn giữa các quan sát. Theo trực giác, thật hợp lý khi mô hình MA có thể được sử dụng để mô tả thành phần bất thường trong chuỗi thời gian về tuổi qua đời của các vị vua Anh, vì chúng ta có thể mong đợi tuổi qua đời của một vị vua cụ thể của Anh có ảnh hưởng nhất định đến tuổi qua đời của một hoặc hai vị vua tiếp theo, nhưng không ảnh hưởng nhiều đến tuổi qua đời của các vị vua sau đó.

5.4.1 Sử dụng hàm `auto.arima()`

Hàm auto.arima() trong gói forecast có thể được sử dụng để tìm mô hình ARIMA thích hợp.

forecast:::auto.arima(kings)

## Series: kings 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -0.7218
## s.e.   0.1208
## 
## sigma^2 = 236.2:  log likelihood = -170.06
## AIC=344.13   AICc=344.44   BIC=347.56

Kết quả cho thấy mô hình ARIMA(0,1,1) là thích hợp để dự báo cho chuỗi tuổi thọ của các vị vua ở Anh dạng gốc hoặc mô hình ARMA(0,1) cho chuỗi sai phân bậc 1.

Một ví dụ khác về bức màn bụi núi lửa ở Bắc bán cầu của Hipel and Mcleod, 1994. Dữ liệu này bao gồm chỉ số màn bụi núi lửa ở bán cầu bắc, từ 1500-1969. Đây là thước đo tác động của bụi và sol khí thải ra môi trường từ các vụ núi lửa phun trào.

volcanodust <- scan("http://robjhyndman.com/tsdldata/annual/dvi.dat", skip = 1)
volcanodustseries <- ts(volcanodust,start = c(1500))
plot.ts(volcanodustseries)

Biểu đồ cho thấy có vẻ như các biến động ngẫu nhiên trong chuỗi thời gian là không đổi theo thời gian, vì vậy một mô hình cộng có thể thích hợp để mô tả biến động của chuỗi này.

Hơn nữa, chuỗi thời gian dường như là dừng ở bậc gốc. Vì level và phương sai của nó dường như gần như không đổi theo thời gian. Vì vậy, chúng ta không cần lấy sai phân chuỗi này để ước lượng với mô hình ARIMA. Ở đây, bậc tích hợp d chính là bằng 0.

Bây giờ chúng ta có thể vẽ biểu đồ tương quan và biểu đồ tương quan từng phần cho 20 độ trễ đầu tiên để lựa chọn một mô hình ARIMA phù hợp.

acf(volcanodustseries, lag.max = 20, na.action = na.pass)

acf(volcanodustseries, lag.max = 20, plot = FALSE)

## 
## Autocorrelations of series 'volcanodustseries', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.666  0.374  0.162  0.046  0.017 -0.007  0.016  0.021  0.006  0.010 
##     11     12     13     14     15     16     17     18     19     20 
##  0.004  0.024  0.075  0.082  0.064  0.039  0.005  0.028  0.108  0.182

Từ biểu đồ tương quan, chúng tôi thấy rằng tự tương quan ở các độ trễ 1, 2 và 3 vượt quá giới hạn ý nghĩa thống kê và tự tương quan giảm dần về 0 sau 3 độ trễ.

Tự tương quan cho độ trễ 1, 2, 3 là dương và giảm độ lớn khi độ trễ tăng (độ trễ 1: 0,666, độ trễ 2: 0,374, độ trễ 3: 0,162).

Tự tương quan cho độ trễ 19 và 20 cũng vượt quá giới hạn ý nghĩa thống kê, nhưng có khả năng điều này là do ngẫu nhiên. Sự tự tương quan ở các độ trễ 4-18 không vượt quá giới hạn đáng kể và chúng ta mong đợi 1 trong 20 độ trễ vượt quá giới hạn ý nghĩa 95% chỉ do ngẫu nhiên.

pacf(volcanodustseries, lag.max = 20, na.action = na.pass)

pacf(volcanodustseries, lag.max = 20, plot = FALSE)

## 
## Partial autocorrelations of series 'volcanodustseries', by lag
## 
##      1      2      3      4      5      6      7      8      9     10     11 
##  0.666 -0.126 -0.064 -0.005  0.040 -0.039  0.058 -0.016 -0.025  0.028 -0.008 
##     12     13     14     15     16     17     18     19     20 
##  0.036  0.082 -0.025 -0.014  0.008 -0.025  0.073  0.131  0.063

Từ biểu đồ tự tương quan từng phần, chúng ta thấy rằng tự tương quan từng phần ở độ trễ 1 là dương và vượt quá giới hạn ý nghĩa thống kê (0,666). Trong khi tự tương quan từng phần ở độ trễ 2 là âm và cũng vượt quá giới hạn ý nghĩa thống kê (-0,126). Sự tự tương quan từng phần giảm dần về 0 sau 2 độ trễ đầu tiên.

Vì biểu đồ tương quan tắt dần về 0 sau độ trễ 3 và biểu đồ tương quan từng phần cũng tắt dần sau 2 độ trễ gợi ý rằng một số mô hình ARMA sau có thể được xem xét: - Mô hình ARMA(2,0), do PACF bằng 0 sau độ trễ 2 và ACF tắt dần về 0 sau 3 độ trễ. - Mô hình ARMA(0,3), vì ACF bằng 0 sau độ trễ 3 và PACF tắt dần về 0. - Mô hình hỗn hợp ARMA(p,q), vì biểu đồ ACF và PACF tắt dần về 0.

5.4.2 Sử dụng tiêu chí thông tin

Chúng ta có thể sử dụng hàm auto.arima() trong gói forecast để tìm mô hình ARIMA phù hợp.

forecast:::auto.arima(volcanodust)

## Series: volcanodust 
## ARIMA(1,0,2) with non-zero mean 
## 
## Coefficients:
##          ar1     ma1     ma2     mean
##       0.4723  0.2694  0.1279  57.5178
## s.e.  0.0936  0.0969  0.0752   8.4883
## 
## sigma^2 = 4897:  log likelihood = -2661.84
## AIC=5333.68   AICc=5333.81   BIC=5354.45

forecast:::auto.arima(volcanodust, ic = "bic")

## Series: volcanodust 
## ARIMA(2,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2     mean
##       0.7533  -0.1268  57.5274
## s.e.  0.0457   0.0458   8.5958
## 
## sigma^2 = 4901:  log likelihood = -2662.54
## AIC=5333.09   AICc=5333.17   BIC=5349.7

Kết quả cho thấy mô hình ARIMA(1,0,2) là mô hình được lựa chọn để mô tả chuỗi thời gian của chỉ số màn bụi núi lửa. Nếu chúng ta chỉ định chọn mô hình theo tiêu chí thông tin BIC thì ARIMA(2,0,0) là mô hình phù hợp.

ARIMA(2,0,0) là một mô hình tự hồi quy bậc hai, hay AR(2). Mô hình này có thể được viết lại dưới dạng tường minh như sau:

\[X_t - \mu = \beta_1 (X_{t-1} - \mu) + \beta_2 (X_{t-2} - \mu) + Z_t\] Trong đó:

$X_t$ là một chuỗi dừng, trong ví dụ là dừng ở bậc gốc
$\mu$ là giá trị trung bình của chuỗi $X_t$,
$Z_t$ là một nhiễu trắng (white noise) với trung bình bằng 0 và phương sai không đổi,
$\beta_1$ và $\beta_2$ là hai tham số cần được ước lượng.

Mô hình tự hồi quy AR thường được sử dụng để lập mô hình chuỗi thời gian cho thấy sự phụ thuộc dài hạn hơn giữa các quan sát liên tiếp. Trong ví dụ này có thể sử dụng một mô hình AR để mô tả chuỗi thời gian của chỉ số màn bụi núi lửa. Như chúng ta kì vọng mức bụi núi lửa và sol khí trong một năm sẽ ảnh hưởng đến chúng trong nhiều năm sau đó, vì bụi và sol khí khó có thể biến mất nhanh chóng.

5.5 Dự báo sử dụng mô hình ARIMA

Khi bạn đã chọn mô hình ARIMA(p,d,q) ứng cử viên tốt nhất cho dữ liệu chuỗi thời gian của mình, bạn có thể ước tính các tham số và sử dụng mô hình đó làm mô hình dự báo để đưa ra dự đoán cho các giá trị tương lai của chuỗi.

Bạn có thể ước tính các tham số của mô hình ARIMA(p,d,q) bằng hàm arima() trong R.

Trở lại ví dụ về tuổi thọ các vị vua ở Anh thì ARIMA(0,1,1) là mô hình được chọn. Để ước lượng mô hình ARIMA này trên R tôi sử dụng hàm arima() như sau:

kingstimeseriesarima <- arima(kingstimeseries, order = c(0,1,1))
kingstimeseriesarima

## 
## Call:
## arima(x = kingstimeseries, order = c(0, 1, 1))
## 
## Coefficients:
##           ma1
##       -0.7218
## s.e.   0.1208
## 
## sigma^2 estimated as 230.4:  log likelihood = -170.06,  aic = 344.13

Ở đây, bậc của các thành phần được thiết lập ở tùy chọn order = c(0,1,1).

Nếu viết cho chuỗi gốc ban đầu thì mô hình ước lượng là ARIMA(0,1,1). Nếu viết cho chuỗi sai phân bậc 1 thì mô hình ước lượng là ARMA(0,1). Một mô hình ARMA(0,1) có thể được viết dưới dạng toán học như sau:

\[X_t - \mu = Z_t - \theta Z_{t-1}\]

Trong đó:

$X_t$ là một chuỗi dừng, trong ví dụ là dừng ở bậc gốc
$\mu$ là giá trị trung bình của chuỗi $X_t$,
$Z_t$ là một nhiễu trắng (white noise) với trung bình bằng 0 và phương sai không đổi,
$\theta$ là tham số cần được ước lượng.

Kết quả cho thấy giá trị ước tính của $\theta$ bằng -0.7218 (kí hiệu là ma1 trong kết quả).

5.6 Tính khoảng tin cậy dự báo

Bạn có thể chỉ định khoảng tin cậy cho khoảng dự báo bằng hàm forecast.Arima() bằng cách thiết lập tùy chọn level của câu lệnh. Ví dụ, để tính khoảng dự báo 99.5%, tôi thiết lập như sau:

forecast:::forecast.Arima(kingstimeseriesarima, h = 5, level = c(99.5))

##    Point Forecast  Lo 99.5  Hi 99.5
## 43       67.75063 25.13941 110.3618
## 44       67.75063 23.52078 111.9805
## 45       67.75063 21.95932 113.5419
## 46       67.75063 20.44939 115.0519
## 47       67.75063 18.98618 116.5151

Chúng ta có thể sử dụng mô hình ARIMA để thực hiện dự báo các giá trị tương lai của chuỗi bằng cách sử dụng hàm forecast.Arima() với thiết lập số kỳ dự báo ở tùy chọn h như sau:

library("forecast")
kingstimeseriesforecasts <- forecast:::forecast.Arima(kingstimeseriesarima, h = 5)
kingstimeseriesforecasts

##    Point Forecast    Lo 80    Hi 80    Lo 95     Hi 95
## 43       67.75063 48.29647 87.20479 37.99806  97.50319
## 44       67.75063 47.55748 87.94377 36.86788  98.63338
## 45       67.75063 46.84460 88.65665 35.77762  99.72363
## 46       67.75063 46.15524 89.34601 34.72333 100.77792
## 47       67.75063 45.48722 90.01404 33.70168 101.79958

Chuỗi thời gian ban đầu của các vị vua Anh bao gồm tuổi qua đời của 42 vị vua Anh. Hàm forecast.Arima() cho ta dự báo tuổi thọ của 5 vị vua tiếp theo, cụ thể là các vị vua từ triều đại 43-47, cũng như tính toán khoảng dự đoán 80% và 95% cho những dự đoán đó.

Tuổi thọ của vị vua thứ 42 của Anh là 56 tuổi, và mô hình ARIMA đưa ra dự báo tuổi thọ của 5 vị vua tiếp theo là 67,8 tuổi. Chúng ta có thể vẽ các tuổi thọ quan sát được của 42 vị vua trong mẫu, cũng như các tuổi thọ dự báo của 5 vị vua tiếp theo sử dụng mô hình ARIMA(0,1,1) như sau:

forecast:::plot.forecast(kingstimeseriesforecasts)

Như trong trường hợp các mô hình làm trơn hàm mũ, nên xem xét liệu các sai số dự báo của mô hình ARIMA có phân phối chuẩn với giá trị trung bình bằng 0 và phương sai không đổi, cũng như còn tồn tại sự tương quan giữa các sai số dự báo hay không.

Ví dụ: chúng ta có thể tạo một biểu đồ ACF của các sai số dự báo từ mô hình ARIMA(0,1,1) cho chuỗi tuổi thọ của các vị vua và thực hiện kiểm định Ljung-Box cho 20 độ trễ như sau:

acf(kingstimeseriesforecasts$residuals, lag.max = 20, na.action = na.pass)

Box.test(kingstimeseriesforecasts$residuals, lag = 20, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  kingstimeseriesforecasts$residuals
## X-squared = 13.584, df = 20, p-value = 0.8509

Vì biểu đồ ACF cho thấy rằng không có sự tự tương quan nào trong 20 độ trễ đầu tiên vượt quá giới hạn ý nghĩa thống kê và giá trị p của kiểm định Ljung-Box là 0,9. Vì vậy, chúng ta có thể kết luận rằng có có rất ít bằng chứng cho sự tự tương quan khác không trong các sai số dự báo ở 20 độ trễ đầu tiên.

Để kiểm tra xem các sai số dự báo có phân phối chuẩn với giá trị trung bình bằng 0 và phương sai không đổi hay không, chúng ta có thể vẽ biểu đồ thời gian và biểu đồ phân phối chuẩn của các sai số dự báo:

plot.ts(kingstimeseriesforecasts$residuals)

plotForecastErrors(kingstimeseriesforecasts$residuals)

Biểu đồ thời gian của các sai số dự báo trong mẫu cho thấy phương sai của các sai số dự báo dường như gần như không đổi theo thời gian (mặc dù có lẽ có phương sai cao hơn một chút đối với nửa sau của chuỗi). Biểu đồ của chuỗi thời gian cho thấy các sai số dự báo có phân phối xấp xỉ phân phối chuẩn với giá trị trung bình gần bằng không. Vì vậy, có thể kết luận các sai số thỏa mãn giả định phân phối chuẩn và phương sai không đổi.

Vì các sai số dự báo liên tiếp không tương quan lẫn nhau và có phân phối xấp xỉ phân phối chuẩn với giá trị trung bình bằng 0 và phương sai không đổi nên mô hình ARIMA(0,1,1) được xem là mô hình dự đoán thích hợp về tuổi thọ của các vị vua Anh.

Tương tự đối với chuỗi về chỉ số bụi núi lửa ở trên. Mô hình ARIMA(2,0,0) là thích hợp để mô tả chuỗi. Để ước lượng mô hình ARIMA(2,0,0) này, tôi sử dụng hàm arima như sau:

volcanodustseriesarima <- arima(volcanodustseries, order = c(2,0,0))
volcanodustseriesarima

## 
## Call:
## arima(x = volcanodustseries, order = c(2, 0, 0))
## 
## Coefficients:
##          ar1      ar2  intercept
##       0.7533  -0.1268    57.5274
## s.e.  0.0457   0.0458     8.5958
## 
## sigma^2 estimated as 4870:  log likelihood = -2662.54,  aic = 5333.09

Như đề cập ở trên, mô hình ARIMA (2,0,0) có thể được viết lại như sau:

\[X_t - \mu = \beta_1 (X_{t-1} - \mu) + \beta_2 (X_{t-2} - \mu) + Z_t\] với $\beta_1$ và $\beta_2$ là hai tham số cần được ước lượng. Kết quả cho thấy giá trị ước lượng của hai tham số này lần lượt là 0.7533 và -0.1268 (kí hiệu là ar1 và ar2 trong kết quả ước lượng của hàm arima).

Bây giờ chúng ta đã ước lượng mô hình ARIMA(2,0,0), do đó chúng ta có thể sử dụng hàm forecast.Arima() để thực hiện dự báo các chỉ số bụi núi lửa trong tương lai. Dữ liệu ban đầu gồm các chỉ số bụi trong giai đoạn 1500-1969. Để thực hiện dự báo chỉ số bụi cho 31 năm sau đó (từ 1970-2000), tôi thực hiện như sau:

volcanodustseriesforecasts <- forecast:::forecast.Arima(volcanodustseriesarima, h = 31)
volcanodustseriesforecasts

##      Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 1970       21.48131 -67.94860 110.9112 -115.2899 158.2526
## 1971       37.66419 -74.30305 149.6314 -133.5749 208.9033
## 1972       47.13261 -71.57070 165.8359 -134.4084 228.6737
## 1973       52.21432 -68.35951 172.7881 -132.1874 236.6161
## 1974       54.84241 -66.22681 175.9116 -130.3170 240.0018
## 1975       56.17814 -65.01872 177.3750 -129.1765 241.5327
## 1976       56.85128 -64.37798 178.0805 -128.5529 242.2554
## 1977       57.18907 -64.04834 178.4265 -128.2276 242.6057
## 1978       57.35822 -63.88124 178.5977 -128.0615 242.7780
## 1979       57.44283 -63.79714 178.6828 -127.9777 242.8634
## 1980       57.48513 -63.75497 178.7252 -127.9356 242.9059
## 1981       57.50627 -63.73386 178.7464 -127.9145 242.9271
## 1982       57.51684 -63.72330 178.7570 -127.9040 242.9376
## 1983       57.52212 -63.71802 178.7623 -127.8987 242.9429
## 1984       57.52476 -63.71538 178.7649 -127.8960 242.9456
## 1985       57.52607 -63.71407 178.7662 -127.8947 242.9469
## 1986       57.52673 -63.71341 178.7669 -127.8941 242.9475
## 1987       57.52706 -63.71308 178.7672 -127.8937 242.9479
## 1988       57.52723 -63.71291 178.7674 -127.8936 242.9480
## 1989       57.52731 -63.71283 178.7674 -127.8935 242.9481
## 1990       57.52735 -63.71279 178.7675 -127.8934 242.9481
## 1991       57.52737 -63.71277 178.7675 -127.8934 242.9482
## 1992       57.52738 -63.71276 178.7675 -127.8934 242.9482
## 1993       57.52739 -63.71275 178.7675 -127.8934 242.9482
## 1994       57.52739 -63.71275 178.7675 -127.8934 242.9482
## 1995       57.52739 -63.71275 178.7675 -127.8934 242.9482
## 1996       57.52739 -63.71275 178.7675 -127.8934 242.9482
## 1997       57.52739 -63.71275 178.7675 -127.8934 242.9482
## 1998       57.52739 -63.71275 178.7675 -127.8934 242.9482
## 1999       57.52739 -63.71275 178.7675 -127.8934 242.9482
## 2000       57.52739 -63.71275 178.7675 -127.8934 242.9482

Chúng ta có thễ vẽ các giá trị của chuỗi ban đầu và các giá trị dự báo như sau:

forecast:::plot.forecast(volcanodustseriesforecasts)

Một điều đáng lo ngại là mô hình đã dự đoán các giá trị âm cho chỉ số bụi núi lửa, nhưng biến này chỉ có thể ghi nhận các giá trị dương! Lý do là các hàm arima() và forecast.Arima() không biết rằng biến chỉ có thể nhận giá trị dương. Rõ ràng, đây không phải là điều mong muốn của mô hình dự báo hiện tại của chúng tôi.

Một lần nữa, chúng ta nên xem xét liệu các sai số dự báo có tương quan với nhau hay không và liệu chúng có phân phối chuẩn với trung bình bằng 0 và phương sai không đổi hay không. Để kiểm tra sự tự tương quan giữa các sai số dự báo liên tiếp, chúng ta có thể tạo biểu đồ tương quan và sử dụng kiểm định Ljung-Box:

acf(volcanodustseriesforecasts$residuals, lag.max = 20, na.action = na.pass)

Box.test(volcanodustseriesforecasts$residuals, lag = 20, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  volcanodustseriesforecasts$residuals
## X-squared = 24.364, df = 20, p-value = 0.2268

Biểu đồ tương quan cho thấy tự tương quan mẫu ở độ trễ 20 vượt quá mức giới hạn ý nghĩa thống kê. Tuy nhiên, điều này có thể là do ngẫu nhiên. Hơn nữa, giá trị p của Kiểm định Ljung-Box là 0,2, chỉ ra rằng có rất ít bằng chứng về tự tương quan khác không trong các sai số dự báo cho 20 độ trễ đầu tiên.

Để kiểm tra xem các sai số dự báo có phân phối chuẩn với giá trị trung bình bằng 0 và phương sai không đổi hay không, chúng ta vẽ biểu đồ thời gian của các sai số dự báo:

plot.ts(volcanodustseriesforecasts$residuals)

plotForecastErrors(volcanodustseriesforecasts$residuals)

Biểu đồ thời gian của các sai số dự báo cho thấy các sai số dự báo dường như có phương sai không đổi theo thời gian. Tuy nhiên, chuỗi thời gian của các sai số dự báo dường như có trung bình âm, chứ không phải là trung bình bằng không. Chúng tôi có thể xác nhận điều này bằng cách tính giá trị trung bình của các sai số là khoảng -0,22.

mean(volcanodustseriesforecasts$residuals)

## [1] -0.2205417

Biểu đồ histogram của các sai số dự báo (ở trên) cho thấy rằng mặc dù giá trị trung bình của sai số dự báo là âm, phân bố của sai số dự báo hơi lệch về bên phải so với đường cong phân phối chuẩn. Do đó, chúng ta không thể thoải mái kết luận rằng sai số dự báo có phân phối xấp xỉ phân phối chuẩn với giá trị trung bình bằng 0 và phương sai không đổi! Như vậy, rất có thể rằng mô hình ARIMA(2,0,0) của chúng ta cho chuỗi thời gian của chỉ số màn bụi núi lửa không phải là mô hình tốt nhất mà chúng tôi có thể xây dựng và gần như chắc chắn mô hình này có thể được cải thiện!

6 Tài liệu tham khảo

Nội dung được phát hành theo điều kiện Creative Commons Attribution-ShareAlike 3.0 Unported.

Phân tích chuỗi thời gian với R - Cơ bản

Vietlod

12 December, 2022