›
›

Phân tích chuỗi thời gian và các kỹ thuật dự báo

pdf 76 trang vanle 2290
Download
Bạn đang xem 20 trang mẫu của tài liệu "Phân tích chuỗi thời gian và các kỹ thuật dự báo", để tải tài liệu gốc về máy bạn click vào nút DOWNLOAD ở trên

Tài liệu đính kèm:

  • pdfphan_tich_chuoi_thoi_gian_va_cac_ky_thuat_du_bao.pdf

Nội dung text: Phân tích chuỗi thời gian và các kỹ thuật dự báo

  1. TRƯỜNG ĐẠI HỌC THỦY LỢI PHÂN TÍCH CHUỖI THỜI GIAN VÀ CÁC KỸ THUẬT DỰ BÁO [Tài liệu giảng dạy ở bậc đại học] Nguyễn Thị Vinh HÀ NỘI 2010 1
  2. MỤC LỤC 1 CHƯƠNG 1: CÁC KHÁI NIỆM CHUNG VỀ DỰ BÁO 1 1.1 Bài toán dự báo 1 1.1.1 Các bài toán 1 1.1.2 Dự báo hỗ trợ quá trình ra quyết định trong các tình huống 1 1.1.3 Tiến trình dự báo chung 2 1.2 Một số khái niệm cơ bản trong dự báo 2 1.2.1 Chuỗi thời gian (Time Series) 2 1.2.2 Các phương pháp hiển thị chuỗi thời gian 3 1.2.3 Các định dạng dữ liệu 4 1.3 Tiêu chuẩn dự báo 6 1.3.1 Các đặc tính thống kê: 6 1.3.2 Các đặc tính định dạng 6 1.4 Liên hệ giữa tính toán hồi qui và dự báo chuỗi thời gian 6 1.5 BÀI TẬP CHƯƠNG 1 7 2 CHƯƠNG 2: CÁC MÔ HÌNH TRƠN 8 2.1 Khái niệm chung về các mô hình trơn 8 2.2 Phương pháp ngây thơ (naive) - phương pháp đơn giản nhất: 8 2.3 Các mô hình trơn không có tính mùa (thời vụ) 9 2.3.1 Mô hình trung bình trượt đơn (Moving Average) 9 2.3.2 Mô hình trung bình trượt với trọng số dạng hàm mũ 9 2.3.3 Các mô hình xu thế 11 2.4 Các mô hình trơn có yếu tố thời vụ (mùa) của Winters 17 2.4.1 Các khái niệm chung 17 2.4.2 Mô hình Winters cho dạng xu thế tuyến tính, thời vụ cộng tính 18 2.4.3 Mô hình Winters cho dạng xu thế mũ, thời vụ nhân tính 18 2.4.4 Mô hình Winters cho dạng xu thế tuyến tính, thời vụ nhân tính (dạng phổ biến nhất) 18 2.4.5 Mô hình Winters cho dạng xu thế mũ, thời vụ cộng tính 19 2.4.6 Các nhận xét chung về các mô hình Winters: 19 2.5 Các phương pháp phân ly (Decomposition) 22 2.5.1 Các công thức chung 22 2.5.2 Phương pháp phân ly cổ điển (Classical Decomposition) 23 2.5.3 Các ví dụ 23 2.6 BÀI TẬP CHƯƠNG 2 26 3 CHƯƠNG 3 : PHÂN TÍCH CHUỖI THỜI GIAN VÀ CÁC MÔ HÌNH CỦA BOX-JENKINS 28 3.1 Các mô hình chuỗi thời gian ARMA (AutoRegressive-Moving Average) 28 3.1.1 Mô hình tự hồi quy bậc p - AR(p) 28 3.1.2 Mô hình trung bình trượt bậc q - MA(q) 29 3.1.3 Mô hình hỗn hợp tự hồi quy-trung bình trượt bậc (p,q) - ARMA(p,q) 29 1
  3. 3.2 Các điều kiện cần về tính dừng và tính khả nghịch 29 3.2.1 Điều kiện dừng 29 3.2.2 Điều kiện khả nghịch 30 3.3 Các trợ giúp cho việc phân tích chuỗi thời gian 31 3.3.1 Biểu diễn đồ họa chuỗi thời gian 31 3.3.2 Hệ số tự tương quan ACF (Auto Correlation Function) 31 3.3.3 Hàm tự tương quan riêng phần PACF 33 3.3.4 Thống kê Q của Box-Pierce 36 3.4 Các ứng dụng của các hệ số tự tương quan 37 3.4.1 Kiểm tra tính ngẫu nhiên của dữ liệu và phần dư 37 3.4.2 Xác định tính dừng của chuỗi thời gian 37 3.4.3 Loại bỏ tính không dừng của chuỗi thời gian 39 3.4.4 Nhận biết tính thời vụ trong chuỗi thời gian 40 3.5 Các mô hình ARIMA 43 3.5.1 Các mô hình ARIMA không có tính thời vụ 43 3.5.2 Các mô hình ARIMA có tính thời vụ 46 3.6 BÀI TẬP CHƯƠNG 3 53 4 CHƯƠNG 4: CÁC PHƯƠNG PHÁP DỰ BÁO CỦA BOX-JENKINS 55 4.1 Các khâu chính trong phương pháp Box-Jenkins 55 4.2 Các nguyên tắc lựa chọn mô hình ARIMA(p,d,q) phù hợp 56 4.3 Các hàm dự báo của các mô hình ARMA(p,q) 58 4.3.1 Một số mô hình ARMA thường gặp: 59 4.3.2 Giới hạn cho phép của các dự báo 60 4.4 Các ví dụ minh họa 60 4.5 BÀI TẬP CHƯƠNG 4 64 5 PHỤ LỤC: GIỚI THIỆU PHẦN MỀM DỰ BÁO SIBYL 65 5.1 Môi trường làm việc của Sibyl 65 5.2 Một số phương pháp dự báo trong Sibyl 66 5.2.1 Các phương pháp trung bình trượt 66 5.2.2 Các phương pháp hồi quy tìm đường cong phù hợp với chuỗi dữ liệu (Trend-Cycle Regression Curve-Fitting Methods) 66 5.2.3 Các phương pháp làm trơn dạng mũ 67 5.2.4 Các phương pháp phân ly 68 5.2.5 Phương pháp Box-Jenkins 69 2
  4. 1 CHƯƠNG 1: CÁC KHÁI NIỆM CHUNG VỀ DỰ BÁO Dự báo là quá trình tạo ra các nhận định về các hiện tượng mà thông thường các đầu ra của chúng còn chưa quan sát được. 1.1 Bài toán dự báo 1.1.1 Các bài toán Dự báo là một trong những yếu tố quan trọng nhất trong việc ra các quyết định quản lý bởi vì ảnh hưởng sau cùng của một quyết định thường phụ thuộc vào sự tác động của các nhân tố không thể nhìn thấy tại thời điểm ra quyết định. Vai trò của dự báo là nhậy cảm trong các lĩnh vực như tài chính, nghiên cứu thị trường, lập kế hoạch sản xuất, hành chính công, điều khiển quá trình sản xuất hay nghiên cứu, Trong giới doanh nhân, các câu hỏi thường xuyên được đưa ra là: Lượng hàng sẽ bán trong tháng tới là bao nhiêu? Tháng này nên đặt mua bao nhiêu hàng? Nên giữ bao nhiêu cổ phiếu ? Nên mua bao nhiêu nguyên liệu ? Mục tiêu bán hàng sắp tới là gì? Có nên tăng nhân công không? 1.1.2 D ự báo hỗ trợ quá trình ra quyết định trong các tình huống i> Điều tiết nguồn tài nguyên sẵn có: Dự báo nhu cầu cho sản phẩm, nguyên liệu, nhân công, tài chính hay dịch vụ như là một đầu vào thiết yếu để điều tiết kế hoạch sản xuất, vận tải, tiền vốn và nhân lực. ii> Yêu cầu thêm tài nguyên: Dự báo giúp xác định tài nguyên cần có trong tương lai (như nhân lực, máy móc thiết bị, vốn ) iii> Thiết kế, lập quy hoạch: Dự báo các hiện tượng thiên nhiên như lũ lụt, hạn hán để thiết kế các công trình như đê, đập, hồ chứa và quy hoạch vùng sản xuất. Nhược điểm của dự báo là không thể tránh khỏi sai số. Trên quan điểm thực tiễn, cần hiểu rõ cả mặt mạnh lẫn mặt hạn chế của các phương pháp dự báo và tính đến chúng trong khi sử dụng dự báo. 1
  5. 1.1.3 Tiến trình dự báo chung Nhận dạng mục đích dự báo Thu thập dữ liệu có liên quan trước thời điểm cần dự báo Biểu diễn đồ hoạ dữ liệu, nhận dạng bất kì dạng mẫu nào Lựa chọn mô hình dự báo phù hợp với dạng dữ liệu và dự báo Tính sai số dự báo cho các giá trị tham số khác nhau và lựachọn tham số thích Áp dụng mô hình đã chọn và phát ra các dự báo cần có Sử dụng các thông tin về chất lượng để chỉnh sửa dự Tốt Chưatốt Sử Đánh giá các sai số dự báo dụng 1.2 Một số khái niệm cơ bản trong dự báo 1.2.1 Chuỗi thời gian (Time Series) Chuỗi thời gian là một dãy dữ liệu được quan sát ở các thời điểm kế tiếp nhau với cùng một đơn vị đo mẫu. Trong chuỗi thời gian, trình tự thời gian đóng một vai trò thực sự quan trọng, vì vậy các tính toán thống kê thông thường như trung bình mẫu, độ lệch quân phương mẫu, khoảng tin cậy, kiểm định các giả thuyết, không còn thích hợp Một chuỗi thời gian thường bao gồm những thành phần sau đây i>. Thành phần ổn định ii>. Thành phần xu thế 2
  6. iii> Thành phần mùa (thời vụ) iv> Thành phần ngẫu nhiên v> Thành phần chu kì (dài hạn) 1.2.2 Các phương pháp hiển thị chuỗi thời gian Phân tích chuỗi thời gian bao gồm việc nghiên cứu dạng dữ liệu trong quá khứ và giải thích các đặc điểm chính của nó. Một trong các phương pháp đơn giản và hiệu quả nhất là hiển thị trực quan chuỗi đó. Các đặc điểm không dễ thấy trong bảng dữ liệu thường nổi lên qua các minh họa đồ thị. t 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 xt 265 275 282 290 292 300 310 318 330 338 347 350 360 365 370 376 382 387 xt/xt-1 104 103 103 101 103 103 103 104 102 103 101 103 101 101 102 102 101 xt-xt-1 10 7 8 2 8 10 8 12 8 9 3 10 5 5 6 6 5 Ba loại đồ thị minh họa chuỗi thời gian là i> Đồ thị của xt theo t: cung cấp lịch sử dữ liệu gốc chưa bị chuyển đổi qua bất cứ phép biến đổi nào, giúp cho việc nghiên cứu xu thế và nhận dạng. 3
  7. ii> Đồ thị của xt/ / xt-1 x 100 theo t: mỗi điểm trên đồ thị này cho biết giá trị hiện thời của chuỗi tăng hay giảm so với giá trị trước đó. Ví dụ giá trị tại thời điểm t = 2 là 102,9% chỉ ra rằng chuỗi đã tăng 2,9% từ thời điểm t = 2 sang thời điểm t = 3. Nếu mọi giá trị đều lớn hơn 100% nhưng theo xu thế giảm dần thì đồ thị đó chứng tỏ rằng chuỗi này có xu thế tăng nhưng tỉ lệ tăng lại giảm dần. (xt/xt-1)% xt / xt-1(%) ~ t 104 103.5 103 102.5 102 101.5 101 100.5 t 0 2 4 6 8 101214161820 iii> Đồ thị của xt – xt-1 theo t: Đồ thị này biểu diễn sự thay đổi giữa các bước thời gian kế tiếp nhau. Nhìn vào đồ thị ta thấy được khoảng các giá trị biến đổi giữa các bước kề nhau. Ví dụ, từ bảng các giá trị xt ~ t ở trang trước, người ta vẽ được 3 đồ thị tương ứng ở các phần i>, ii>, iii>. 1.2.3 Các định dạng dữ liệu Trước khi áp dụng bất cứ một phương pháp dự báo khoa học cho một tình huống nào, cần phải ghép nối các thông tin (dữ liệu có liên quan) về tình huống đó càng nhiều càng tốt. Những dữ liệu đó được phân thành 2 loại: i> Các dữ liệu bên trong, ví dụ số liệu sản phẩm bán ra trong quá khứ, ii> Các dữ liệu bên ngoài, ví dụ như các thống kê của ngân hàng về tình hình tài chính của công ty (phản ánh thông tin bên trong). 4
  8. Từ các thông tin này, người làm dự báo phải chọn ra thông tin liên quan nhiều nhất đến tình huống cần dự báo. Chẳng hạn, trong dự báo bán hàng, báo cáo hàng bán được trong quá khứ của công ty sẽ cung cấp những thông tin tối thiểu cho việc dự báo. Thông tin tối thiểu cần thỏa mãn các yêu cầu về: - Tính liên quan: Nó có phải là thông tin liên quan trực tiếp nhất không? - Độ tin cậy: Dữ liệu được thu thập như thế nào? Có đáng tin cậy không? - Tính thời sự: Liệu các thông tin mới nhất đã được cập nhật chưa? Chúng có sẵn khi cần không? Thời vụ Ổn định (không có xu thế ) (trung bình và phương sai không đổi) Xu thế tuyến tính tăng Xu thế tuyến tính giảm Chu kì dài hạn Xu thế tuyến tính tăng và thời vụ nhân tính Xu thế tuyến tính tăng và thời vụ cộng tính Một số định dạng dữ liệu Khi đã có những thông tin tối thiểu cần thiết, ta cần phải nghiên cứu đặc điểm của nó bằng cách minh họa đồ thị. Dạng dữ liệu quá khứ là rất quan trọng vì nó quyết định 5
  9. việc lựa chọn mô hình dự báo. Mô hình dự báo được chọn phải tương thích với dạng dữ liệu mẫu trong quá khứ. 1.3 Tiêu chuẩn dự báo Các tiêu chuẩn chung đánh giá sự thành công của một mô hình dự báo khi áp dụng vào một tập dữ liệu là: i> Trùng càng nhiều với các thay đổi ngẫu nhiên trong dữ liệu càng tốt. ii> Không vượt quá xa bất kì một đặc tính nào của dữ liệu Xét về mặt sai số, hai loại đặc tính cần quan tâm khi thử nghiệm một công thức dự báo trên dữ liệu là 1.3.1 Các đặc tính thống kê: Một phương pháp dự báo tốt thường cho sai số trung bình nhỏ. Trong các mô hình dự báo, người ta thường sử dụng các loại sai số như 1 MAE = e (Mean Absolute Error) n ∑ i 1 MSE = e2 (Mean Square Error) n ∑ i RMSE = MSE (squareRoot Mean Square Error) ở đây sai số ei =xi – fi với fi là dự báo của xi 1.3.2 Các đặc tính định dạng Trong các mô hình dự báo, sự có mặt của các dạng sai số (như tính lệch, tính chu kì, tính kiên định, ) đều bị xem là dấu hiệu không tốt. Sự xuất hiện của bất cứ xu thế nào trong sai số cũng nên khử càng nhanh càng tốt. Có thể sai phân hóa chuỗi các giá trị ban đầu để đối phó với các tác động này Tóm lại có hai tiêu chuẩn dự báo về định lượng và định tính là: sai số nhỏ và không tuân theo một định dạng nào. 1.4 Liên hệ giữa tính toán hồi qui và dự báo chuỗi thời gian Tính toán hồi qui dựa trên quan hệ nhân – quả của hệ thống và cực tiểu sai số bằng phương pháp bình phương bé nhất Dự báo chuỗi thời gian dựa trên quan hệ nội tại của dữ liệu để phát ra các dự báo cho các bước thời gian tiếp theo. 6
  10. 1.5 BÀIH TẬP C ƯƠNG 1 1. Trong các định dạng có thể có của chuỗi thời gian, những định dạng nào có tính loại trừ nhau? 2. Giải thích tại sao một kết quả dự báo có sai số không ngẫu nhiên, tức là tuân theo một định dạng nào đó, là một dự báo không tốt? 7
  11. 2 CHƯƠNG 2: CÁC MÔ HÌNH TRƠN 2.1 Khái niệm chung về các mô hình trơn Cơ sở của các phương pháp này là làm trơn (lấy trung bình hoặc trung bình có trọng số) các quan sát trong quá khứ của chuỗi thời gian để nhận được dự báo cho tương lai. Trong việc làm trơn các giá trị quá khứ, các sai số ngẫu nhiên được tính trung bình. Các mô hình trơn dùng trong dự báo thích hợp cho một số tình huống. Các ưu điểm chính của các phương pháp làm trơn là: i> Chi phí thấp ii> Dễ dùng (ở những nơi có thể áp dụng được) iii> Tốc độ tính nhanh (ở những nơi chấp nhận được) Những phương pháp làm trơn rất hấp dẫn khi cần phải dự báo ở rất nhiều bước thời gian tương lai, chẳng hạn trong công tác kiểm kê. 2.2 Phương pháp ngây thơ (naive) - phương pháp đơn giản nhất: Giả sử người quản lý siêu thị muốn biết một khách hàng điển hình tiêu bao nhiêu tiền cho một lần mua sắm. Lấy ngẫu nhiên một mẫu 12 khách hàng và nhận được kết quả sau: Khách hàng 1 2 3 4 5 6 7 8 9 10 11 12 Số tiền đã tiêu ($) 19 18 19 22 21 17 23 19 19 22 21 20 Có thể lấy giá trị làm cực tiểu sai số MSE, trong trường hợp này là trung bình mẫu X = $20. ¾ Các ưu điểm của phương pháp trung bình i> Cực tiểu sai số ii> Ước lượng không chệch iii> Cho dự báo tốt nếu dữ liệu có tính ổn định (trung bình không đổi) và tính ngẫu nhiên (không có xu thế tăng /giảm, không có tính thời vụ hay chu kì). ¾ Các nhược điểm của phương trung bình i> Cho kết quả tồi nếu dữ liệu có tính xu thế hoặc có định dạng xác định ii> Cần mẫu có dung lượng lớn iii>Dự báo tồi nếu có đột biến Kết luận: chỉ sử dụng phương pháp ngây thơ khi chuỗi thời gian có tính ổn định, ngẫu nhiên, và khi không biết phương pháp dự báo nào khác. 8
  12. 2.3 Các mô hình trơn không có tính mùa (thời vụ) 2.3.1 Mô hình trung bình trượt đơn (Moving Average) ¾ Phương pháp: Lấy trung bình N giá trị liên tiếp của các quan sát gần nhất làm dự báo cho thời điểm thứ N+1. Thuật ngữ trung bình trượt có nghĩa là quan sát cũ nhất sẽ bị loại đi mỗi khi có quan sát mới. Nói cách khác, số quan sát trong khi tính là không đổi và chỉ bao gồm các quan sát gần với hiện tại nhất. ¾ Lập công thức: ft+1 = (xt + xt-1 + + xt-N+1) / N (2.1) = (xt-1 + xt-2 + + xt-N) / N + xt/ N – xt-N / N hay ft+1 = ft + xt/ N – xt-N / N (2.2) ¾ Nhận xét: i> Dự báo ở thời điểm t+1 chỉ là điều chỉnh của dự báo ở thời điểm t trước đó. Khi N tăng đủ lớn thì lượng điều chỉnh xt / N – xt-N / N → 0 và trung bình trượt trở thành trung bình mẫu như phương pháp ngây thơ, độ chính xác thấp. ii> Chỉ nên áp dụng phương pháp này khi số giá trị quan sát được là ít và tập dữ liệu có tính ổn định theo thời gian ¾ Ví dụ: Bảng dưới đây cho biết lượng hàng bán ra của các tháng 1, 2, , 11. Nếu sử dụng mô hình trung bình trượt MA với N = 1 ta coi lượng hàng bán ra của tháng trước là dự báo cho tháng sau; với N = 11 ta sử dụng trung bình mẫu cho dự báo của tháng 12; với N = 3 ta sử dụng trung bình của 3 tháng gần nhất làm dự báo cho tháng tới Lương hàng Dự báo bán ra 3 tháng Bán ra 200 Dự báo 135 350 195 300 197.5 176.7 310 175.8 250 175 234.2 155 227.5 200 130 213.3 220 153.3 150 277 168.3 235 209.0 100 -1135791113 - 244.0 2.3.2 Mô hình trung bình trượt với trọng số dạng hàm mũ (Exponentially Weighted Moving Averages) hay mô hình trơn dạng mũ đơn ¾ Phương pháp: Hai hạn chế của mô hình MA là: 9
  13. i> N giá trị quá khứ bắt buộc phải có đủ ii> Trọng số trung bình cho các quan sát là như nhau (1 / N) Trên thực tế, các quan sát càng gần càng chứa nhiều thông tin cho các giá trị sắp xảy ra, do đó cần cho chúng các trọng số lớn hơn so với các quan sát ở xa ¾ Lập công thức: Giả sử chuỗi dữ liệu quan sát được là ổn định (có trung bình không đổi) và không có quan sát thứ N-t . Khi đó từ công thức (2,2) lấy ft thay cho xN-t ta được ft+1 = ft + xt / N – ft / N = (1-1/N) ft + xt / N, vì N > 0 nên 0 Từ công thức dự báo (2.3) → ft+1 = w xt +(1 – w) ft = ft + w (xt – ft ) hay ft+1 = ft + w et (dự báo mới bằng tổng của dự báo cũ và điều chỉnh sai số). Đây chính là nguyên tắc phản hồi hay phương pháp thích ứng của dự báo. ii> Một số trường hợp riêng: w = 0: ft+1 = ft w = 1: ft+1 = xt w ≈ 1: cho các dự báo phản ánh các thay đổi gần đây nhất w = 0,1: ft+1 = 0,1xt + 0,09xt-1 + 0,081xt-2 + cho các dự báo xấp xỉ nhau w = 0,9: ft+1 = 0,9xt + 0,09xt-1 + 0,009xt-2 + dự báo bám theo mẫu 1 bước iii> Chú ý rằng việc phản hồi sự biến đổi của mẫu được cải thiện khi w gần 1. Tuy nhiên việc phản hồi được thực hiện nhanh hay chậm còn tùy vào khả năng làm trơn các dao động ngẫu nhiên. iv> Các ưu điểm của phương pháp EWMA là không cần biết nhiều số liệu quá khứ và tính toán đơn giản; dữ liệu càng gần càng có trọng số lớn; thích hợp khi phải dự báo cho nhiều bước thời gian (khi đó w thường là 0,2 hoặc 0,3) ¾ Một số vấn đề nảy sinh và cách khắc phục: i> Thời điểm đầu tiên t =1: không có dự báo cho thời điểm trước đó để tính dự báo f1 theo công thức (2.3). Các giải pháp là lấy f1 = x1 hoặc f0 = x hoặc sử dụng trung bình cộng của vài giá trị đầu làm giá trị f0; ii> Chọn giá trị w theo một trong ba tiêu chí sau 10
  14. ● w là tốt nhất cho mô hình theo nghĩa sai số MSE là nhỏ nhất. Giá trị này phải được tính thử cho các giá trị w khác nhau để lựa chọn. Trong ví dụ trên w = 0,1 MSE = 3438,3 w = 0,5 MSE = 4347,2 w = 0,9 MSE = 5039,4 Trong ví dụ này, MSE giảm khi w giảm, chứng tỏ dữ liệu là ngẫu nhiên. ● Ở một số bước đầu nên chọn w gần 1 vì không có f0 để tính toán. Có thể tiến hành chọn các wt lớn hơn giá trị tối ưu. Ví dụ khi w = 0,2 là tối ưu thì nên chon wt = 1/t cho đến khi wt < 0,2 . Vậy w1 = 1,0 w2 = 0,5 w3 = 0,33 w4 = 0,25 w = 0,2 với t ≥ 5 ● Định ra các giá trị dường như là tốt nhất cho mỗi tình huống cụ thể, chẳng hạn w=1 khi t = 1; w = 0,3 khi t = 2, 3, 4 và w = 0,3 khi t ≥ 5 2 135 200 200 200 3 195 193.5 167.5 141.5 4 197.5 193.7 181.3 189.7 5 310 194.0 189.4 196.7 6 175 205.6 249.7 298.7 7 155 202.6 212.3 187.4 8 130 197.8 183.7 158.2 9 220 191.0 156.8 132.8 10 277 193.9 188.4 211.3 11 235 202.2 232.7 270.4 12 - 205.5 233.9 238.5 Bán ra w = 0,1 w = 0,5 w = 0 350 300 250 200 150 100 -1 1 3 5 7 9 11 13 2.3.3 Các mô hình xu thế ¾ Đặt vấn đề: Việc áp dụng các mô hình trung bình trượt cho tập dữ liệu chứa xu thế (tăng hoặc giảm) sẽ cho những dự báo thiên nhỏ hoặc thiên lớn so với giá trị thực. Giả sử có N quan sát xt, t = 1, , N theo xu thế tăng tuyến tính như hình vẽ. Ta gọi mức tăng của mẫu tại thời điểm t là 11
  15. mt = a + bt Xu thế tăng tuyến tính trong đó m t = a + bt a = mức tăng tại t = 0 a b = độ dốc Các dự báo được tạo ra tại gốc t = N sẽ là fN+τ = a + b(N + τ), τ = 1, 2, hay fN+τ = m N + b τ Vai trò của các mô hình xu thế là ước lượng mN và b từ các dữ liệu quá khứ. * * Kí hiệu các ước lượng đó là m N và b ta có fN+τ = mN + b τ (2.4) Các cách ước lượng khác nhau cho ta các mô hình tuyến tính khác nhau ¾ Mô hình bình phương bé nhất (Least Mean Square) Việc cực tiểu bình phương tổng các sai số N 2 S=∑ [xt − (a + bt) ] t= 1 dẫn đến các giá trị a* = x − b* t N N N N∑∑∑ txt − t x t (2.5) * t=== 1 t 1 t 1 b = N N N∑∑ t2 − ( t)2 t== 1 t 1 ta có công thức (2.4) với mN = a + b N cho bởi (2.5) ¾ Ví dụ: Cho chuỗi quan sát t 1 2 3 4 5 6 7 8 9 10 xt 60 70 85 60 99 68 106 75 66 124 → f10+τ = 102.81 + 4.58 τ → f11 = 107,4; f12 = 112 Nhận xét: Công thức (2.5) sử dụng trọng số bình quân để tính mN và b ¾ Mô hình trung bình trượt kép (DMA) Lập công thức Mô hình này là sự mở rộng của mô hình MA bằng cách sử dụng các số hạng bám theo xu thế của mẫu. Hai giá trị trung bình trượt được tính tại thời điểm T là 12
  16. ⎧ x+ x + + x 1 N− 1 M = T T− 1 T− N + 1 = x ⎪ T N N ∑ T− i ⎪ i= 0 ⎨ N− 1 ⎪ (2) MT+ M T− 1 + + MT− N + 1 1 MT = = ∑ MT− i ⎩⎪ N N i= 0 N là số bước thời gian được chọn để lấy trung bình trượt MT Ta có công thức tương đương ⎧ x− x MM= + TTN− ⎪ T T− 1 N ⎨ MM− ⎪MM(2) =(2) + TTN− ⎩⎪ T T− 1 N Công thức tính các dự báo tại thời điểm t = n cho τ bước phía trước: * fN+τ = mN + b* τ với ⎧m* = 2M − M(2) , ⎪ N NN (2.6) ⎨ 2 (2) ⎪b* = (MNN - M ) ⎩ N− 1 Chứng minh công thức (2.6) Ở thời điểm T, mức tăng mT = a + bT Ở thời điểm T-1, mức tăng mT-1 = a + b(T–1) = mT – b Ở thời điểm T-2, mức tăng mT-2 = a + b(T–2) = mT - 2b Vậy kì vọng của MT là N− 1 N− 1 N− 1 N− 1 1 1 1 1 E(MT ) = ∑ E(xT− i ) = ∑ mT− i = ∑∑(mT − ib) =(NmT − b i) N i= 0 N i= 0 N i= 0 N i= 0 1 N(N− 1) N− 1 =(Nm − b )= m − b (2.7) N T 2 T 2 1 → E(M(2) ) = E(M+ M + + M ) T N T T− 1 T− N + 1 ⎡ N− 1 N− 1 N− 1 ⎤ =mT − b+ mT− 1 − b+ + mT− N + 1 − b ⎣⎢ 2 2 2 ⎦⎥ 1 N− 1 =[]m + m + + m − b N T T− 1 T− N + 1 2 1 N− 1 =[]m + (m − b) + (m − 2b) + + (m − (N − 1)b) − b N TT T T 2 N− 1 N− 1 =m − b − b= m − (N − 1)b (2.8) T 2 2 T 13
  17. Sử dụng phương pháp ước lượng các moment ta nhận được ⎧ * N− 1 ⎪MTT= m − b * ⎨ 2 ⎪ (2) * ⎩MT = mT − (N − 1)b * Giải hệ 2 phương trình đại số tuyến tính 2 * (2) t xt MT MT ẩn mT và b* ta nhận được (2.6) Ví dụ: Dùng trung bình trượt kép với 1 60 N = 6 tính các dự báo với τ = 1 và τ = 2 2 70 * (2) 3 85 m14 = 2M14 − M14 4 60 = 209,34 – 94,89 = 114,45 5 88 2 b* =(M − M(2) ) 6 66 71.50 5 14 14 7 106 79.17 = 2(104,67-94,89) / 5 = 3,91 8 75 80.00 Vậy f14+τ = 114,45 + 3,91 τ 9 86 80.17 τ = 1 f15 = 118,36; τ = 2 f16= 122,27 10 124 90.83 ¾ Mô hình trơn dạng mũ kép 11 122 96.50 83.03 (DEWMA) 12 87 100.00 87.78 Hạn chế của các mô hình trung bình trượt đơn hay kép là 13 89 97.17 90.78 i> Đòi hỏi N dữ liệu cuối 14 120 104.67 94.89 ii> Trọng số như nhau ở N điểm x 140 này, trọng số 0 cho các điểm khác 120 Phương pháp làm trơn dạng mũ kép sẽ 100 80 khắc phục được các hạn chế trên và trong 60 đa số các trường hợp là thích hợp hơn 40 trung bình trượt kép 20 0 t 051015 Công thức: Gọi xi là dữ liệu gốc ở thời điểm thứ i Si là giá trị làm trơn dạng mũ đơn ở thời điểm thứ i Si’ là giá trị làm trơn dạng mũ kép ở thời điểm thứ i ai là ước lượng của a ở thời điểm thứ i bi là ước lượng của b ở thời điểm thứ i Ta có các quan hệ giữa chúng Si = αXi + (1 – α) Si-1 (2.9) Si’ = αSi + (1 – α) S’i-1 (2.10) 14
  18. Từ đó người ta suy ra được ai = 2Si – Si’ (theo công thức trung bình trượt kép) (2.11) bi = α(Si – Si’) / (1 – α) (2.12) và công thức dự báo DEWMA là fN+τ = aN + bN τ (2.13) Tham số trơn α : Về mặt lí thuyết, α có thể nhận bất cứ giá trị nào giữa 0 và 1. Thực nghiệm cho thấy rằng giá trị tối ưu của α nằm giữa 0,1 và 0,2. α= 0,1cho các dự báo bảo thủ α= 0,2 cho các dự báo phản hồi hệ thống tốt hơn. ' Các giá trị ban đầu a0, b0, S0 và S0 : b0 = x2 – x1 và a0 = x1 – b0 = 2x1 – x2 ● Có thể lấy trung bình N quan sát sau cùng làm ước lượng của S0 x+ x+ + x 1− α S = 1 2 N → SS' = − b 0 N 0 0 α 0 ● Có thể sử dụng các ước lượng thống kê cho a0, b0: chẳng hạn để sử dụng phương pháp làm trơn dạng mũ kép từ chuỗi 11 quan sát, ta có thể lấy hồi qui tuyến tính các giá trị này làm ước lượng mức tăng và độ dốc a0, b0 Khuyến nghị: Phương pháp này thích hợp cho dữ liệu không có yếu tố mùa và không ổn định (có xu thế tăng hoặc giảm) Ví dụ: Cho chuỗi 24 số liệu một mặt hàng bán ra của 24 tháng. Hãy dự báo mức bán ra của tháng tiếp theo với tham số trơn α = 0,2 Bước 1: Sử dụng phương pháp hồi quy tuyến tính cho các dữ liệu quan sát được ta tính được mức tăng và độ dốc cho xu thế chung của mô hình mt = 275 + 10,88 t, t = 1, 2, , 24 → chọn a0 = 275 và b0 = 10,88 x+ x+ + x S = 1 2 24 = 411 0 24 ' 1− α 8,0 → SS0 =0 − b0 = 411 − 10,88= 367,5 α 2,0 t xt St S't at bt et α = 0.2 0 411 367.5 275 10.88 1 317 392 372.4 412 4.944 -100 2 194 353 368.5 337 -3.97 -139 3 312 344 363.7 325 -4.8 -8.4 4 316 339 358.7 319 -4.98 2.13 15
  19. 5 322 335 354 317 -4.65 9.86 6 334 335 350.2 320 -3.78 17.8 7 317 332 346.5 317 -3.75 4.24 8 356 336 344.5 328 -2.02 29.7 9 428 355 346.5 363 2.049 63 10 411 366 350.4 382 3.891 25.6 11 494 392 358.6 425 8.233 61.3 12 412 396 366.1 425 7.404 -21 13 460 409 374.5 443 8.496 8.99 14 395 406 380.8 431 6.256 -42 15 392 403 385.3 421 4.452 -33 16 447 412 390.6 433 5.319 8.56 17 452 420 396.4 443 5.861 2.82 18 571 450 407.2 493 10.73 67.2 19 517 463 418.4 509 11.26 -2.8 20 397 450 424.8 476 6.351 -85 21 410 442 428.3 456 3.473 -50 22 579 470 436.5 503 8.253 68.2 23 473 470 443.2 497 6.741 -31 24 558 488 452.2 523 8.905 25.7 532. f25 = 29 Bước 2: Tính các Si và Si’ theo công thức (2.9) và (2.10), i = 1, 2, , 24 rồi áp dụng công thức (2.11) và (2.12) ta tính được ai, bi và dự báo được f24+1 = a24 + b24 . (1) = 523,4 + 8,9 = 532,3 ≈ 532 Nhận xét: Sai số là đại lượng ngẫu nhiên. ¾ Mô hình Holt Mô hình Holt tương tự như mô hình trơn dạng mũ kép ngoại trừ việc nó không áp dụng công thức trơn kép mà tách riêng việc làm trơn các giá trị xu thế. Điều này làm tăng tính mềm dẻo, vì nó cho phép phần xu thế được làm trơn với tham số khác tham số được sử dụng trong chuỗi quan sát ban đầu. Cụ thể là: ai = αxi + (1 – α) (ai-1 + bi-1) là mức tăng ở thời điểm i bi = β(ai – ai-1) + (1 – β) bi-1 là xu thế (gradient) ở thời điểm i Công thức dự báo: fn+τ = an + bnτ (2.10) Các giá trị ban đầu của a và b là a0 = 2x1 – x2; b0 = x2 –x1 Các giá trị của α, β: Nếu có sẵn một tập các giá trị ban đầu của dữ liệu thì nên sử dụng nó để tìm ra các giá trị α, β tốt nhất. Nếu ta lấy sai số trung bình bình phương (MSE) làm tiêu chuẩn ước lượng, ta có thể ước lượng một khoảng các giá trị khác nhau của α, β. Ví dụ: Cho chuỗi dữ liệu hàng bán ra của 12 tháng năm ngoái. Hãy dự báo mức bán ra của tháng Giêng năm nay với α = 0,2 và β = 0,3 Nhận xét: Nếu số lượng quan sát ít thì các phương pháp dự báo đều cho kết quả nghèo nàn, vì vậy các dự báo nhận được qua vài quan sát ban đầu nên bỏ qua khi tính sai số 16
  20. MSE. Các dạng mô hình trơn bậc cao hơn có thể sử dụng khi xu thế của mẫu có dạng bậc hai, dạng mũ, 0 440 -123 1 317 317 -123 317.00 2 194 194 -123 194.00 quan sát 520 3 312 119 -109 71.00 470 4 316 71.7 -90.2 10.66 420 5 322 49.6 -69.8 -18.49 370 6 334 50.7 -48.5 -20.18 320 7 317 65.1 -29.6 2.11 270 8 356 99.6 -10.4 35.45 220 9 428 157 9.92 89.14 170 120 10 411 216 24.6 166.83 02 468101214 11 494 291 39.8 240.24 12 412 347 44.7 330.78 α =0.2 13 Dư báo 391.7 β =0.3 2.4 Các mô hình trơn có yếu tố thời vụ (mùa) của Winters 2.4.1 Các khái niệm chung Các mô hình này có dạng trơn bậc cao hơn, ưu điểm nổi trội của chúng là sự kết hợp chặt chẽ giữa tính xu thế và yếu tố thời vụ. Các bước phân tích chung Bước 1: Vẽ đồ thị biểu diễn chuỗi thời gian xt ~ t Bước 2: Phân tích ban đầu a) Dữ liệu có thể hiện i> Yêú tố thời vụ? ii> Tính xu thế? b) Nếu có xu thế thì đó là xu thế tuyến tính hay xu thế mũ, có tắt dần không? c) Nếu có yêú tố thời vụ thì đó là tác động cộng tính hay nhân tính, với bước thời vụ là bao nhiêu? Việc nhận dạng dữ liệu sẽ dẫn đến sự lựa chọn mô hình dự báo phù hợp Xu thế tuyến tính tăng và thời vụ nhân tính Xu thế tuyến tính tăng và thời vụ cộng tính xu thế chung xu thế chung 17
  21. Đối với dạng thời vụ cộng tính, các biến đổi theo thời vụ (loại trừ các nhiểu động) là không đổi về mức độ trung bình hoặc xu thế. Đối với dạng này, yếu tố thời vụ thường được ước lượng bởi sự khác biệt giữa giá trị quan sát được với xu thế chung. Ta có mô hình Xt = Tt + It + at trong đó Tt là giá trị xu thế tại thời điểm t, It là yếu tố thời vụ tại thời điểm t và at là nhiễu động tại thời điểm t. Đối với dạng thời vụ nhân tính, do tác động thời vụ tăng / giảm so với mức độ trung bình nên yếu tố thời vụ thường được ước lượng bởi tỉ lệ tăng / giảm so với xu thế chung. Ta có mô hình Xt = Tt It + at hoặc Xt = Tt It at Các mô hình Winters dưới đây đều bao gồm các phương trình trơn dạng mũ tách biệt cho phần xu thế và phần thời vụ 2.4.2 Mô hình Winters cho dạng xu thế tuyến tính, thời vụ cộng tính Các phương trình tính toán các thành phần bao gồm: St = α (xt – It-L) + (1 – α) (St-1 + bt-1) bt = β (St – St-1) + (1 – β) bt-1 (tương tự như mô hình Holt) (2.11) It = γ (Xt – St) + (1 – γ) It-L trong đó St là mức trơn tại thời điểm t bt là xu thế tại thời điểm t It là yếu tố thời vụ tại thời điểm t L là độ dài của thời vụ Dự báo tại thời điểm t = n cho các bước tiếp theo τ = 1, 2, 3, L là fn+τ = Sn + bn τ + In+τ-L 2.4.3 Mô hình Winters cho dạng xu thế mũ, thời vụ nhân tính Các phương trình tính toán các thành phần bao gồm: x t S t =α +(1 − α ) St− 1 b t − 1 I t− L S b t =β +(1 − β ) b t− 1 (2.12) S t− 1 x t I t =γ +(1 − γ ) I t− L S t Dự báo tại thời điểm t = n cho các bước tiếp theo τ = 1, 2, 3, là fn+τ = Sn bn τ In+τ-L 2.4.4 Mô hình Winters cho dạng xu thế tuyến tính, thời vụ nhân tính (dạng phổ biến nhất) Các phương trình tính toán các thành phần bao gồm: 18
  22. x t S t =α +(1 −α ) (St− 1 + b t− 1 ) I t− L (2.13) b t =β (St − S t-1 ) + (1 − β ) b t− 1 x t I t =γ +(1 − γ ) I t− L S t Dự báo tại thời điểm t = n cho các bước tiếp theo τ = 1, 2, 3, là fn+τ = (Sn + bn τ) In+τ-L 2.4.5 Mô hình Winters cho dạng xu thế mũ, thời vụ cộng tính Các phương trình tính toán các thành phần bao gồm: St =α (xt − I t-L ) + (1 − α ) St− 1 b t − 1 S (2.14) bt =β +(1 − β ) bt− 1 St− 1 It=γ (x t − S t ) + (1 − γ ) It− L Dự báo tại thời điểm t = n cho các bước tiếp theo τ = 1, 2, 3, là fn+τ = Sn bn τ + In+τ-L 2.4.6 Các nhận xét chung về các mô hình Winters: ¾ Ưu điểm: Dễ hiểu, sử dụng nhiều trong thực tế, rất phù hợp cho dạng dữ liệu có tính xu thế và yếu tố thời vụ biến đổi. ¾ Nhược điểm : Đòi hỏi 3 tham số trơn, một khi đã được tính toán tối ưu về sai số thì khó điều chỉnh khi nhập thêm quan sát mới. Chú ý : Để tính toán tối ưu các tham số α, β, γ cần tính các giá trị ban đầu S0, b0, và I1, I2, , IL có một số cách sau đây: Cách 1: Dự báo lùi : dùng chuỗi xt. xt-1, , x1 dự báo các giá trị quá khứ x0, x-1, phục vụ cho việc ước lượng S0, b0, và I1, I2, , IL Cách 2: Tách dữ liệu làm 2 phần : ● Phần 1: dùng để ước lượng S0, b0, và I1, I2, , IL. Giả sử có các quan sát cho m thời vụ đầu và x j là trị trung bình của các quan sát ở thời vụ thứ j, với j = 1, 2, , m. Ta có các ước lượng xm− x 1 b = ; 0 ()m− 1 L L S= x1 − b ; 0 2 0 Yếu tố thời vụ tại các thời điểm t = 1, 2, , mL được tính theo công thức x t It = xi − [(L + 1)/2 − j] b0 19
  23. với xi là trị trung bình của thời vụ thứ i tương ứng với thời điểm t, j là vị trí của thời điểm t trong thời vụ thứ i (ví dụ với L+1 ≤ t ≤ 2L thì i = 2, nếu t = L+1 thì j = 1). Lấy trung bình các It trong m thời vụ ta được L giá trị 1 m− 1 It = ∑It+ kL ∀ t = 1, 2, ,L m k= 0 Cuối cùng, các giá trị ban đầu I1, I2, , IL được chọn là chuẩn hóa của các đại lượng It tương ứng It It = L ∀ t = 1, 2, , L 1/L∑ Ik k= 1 ● Phần 2 : dùng để tối ưu hóa α, β, γ theo các mục tiêu làm cực tiểu MSE, RMSE hay MAE. Các kỹ thuật dò tìm có thể là phương pháp thử sai, phương pháp đường dốc nhất, ¾ Ví dụ: Cho dãy 48 số liệu một loại nước giải khát đóng chai bán ra hàng tháng (tính theo kiện) của một hãng trong 4 năm liền. Với các tham số trơn α = 0,2 β = 0,1 và γ = 0,1 hãy sử dụng bảng tính Excel dự báo lượng hàng sẽ bán trong 4 tháng tới. Giải: 1. Đồ thị biểu diễn lượng hàng bán ra theo tháng cho thấy biên độ thời vụ (L=12) tăng theo lượng hàng bình quân bán ra (có xu thế tăng tuyến tính), do đó mô hình Winters với xu thế tuyến tính, thời vụ nhân tính là lựa chọn phù hợp. 2. Số liệu của 2 năm đầu được dùng để tính các giá trị ban đầu, ta có x1 = 359, 42;x2 = 493, 58; 493,58 − 352,75 b = = 12,01 0 (2− 1).12 12 S= 352,75 − .10,49= 289,83 0 2 Cuối cùng các giá trị ban đầu I1, I2, , IL được chọn là chuẩn hóa của các đại lượng It tương ứng It It = L ∀ t = 1, 2, ,L 1/L∑ Ik k= 1 3. Phần còn lại dùng để tối ưu hóa α, β, γ theo các mục tiêu làm cực tiểu MSE, RMSE hay MAE. Các kỹ thuật dò tìm có thể là phương pháp thử sai, phương pháp đường dốc nhất, 20
  24. t xt St bt It ft et 1 143 300.31 10.49 0.48 143.02 -0.02 L = 12 x1TB = 352.75 2 138 293.46 8.75 0.60 191.39 -53.39 α = 0.2 x2TB = 478.58 3 195 301.92 8.72 0.65 195.92 -0.92 β = 0.1 b0 = 10.49 4 225 314.52 9.11 0.69 211.81 13.19 γ = 0.1 S0 = 289.83 5 175 320.06 8.75 0.57 185.21 -10.21 m = 2 6 389 329.78 8.85 1.17 383.32 5.68 7 454 337.80 8.77 1.36 459.66 -5.66 t xt It It TB It ban đầu 8 618 349.58 9.07 1.71 592.22 25.78 1 143 0.48 0.47 0.48 9 770 362.16 9.42 2.05 734.09 35.91 2 138 0.45 0.60 0.62 10 564 388.56 11.12 1.26 459.12 104.88 3 195 0.62 0.64 0.65 11 327 391.31 10.28 0.91 365.20 -38.20 4 225 0.69 0.67 0.68 12 235 402.68 10.39 0.58 231.88 3.12 5 175 0.52 0.56 0.57 13 189 409.83 10.07 0.47 196.71 -7.71 6 389 1.12 1.14 1.17 14 326 444.35 12.51 0.61 252.46 73.54 7 454 1.27 1.33 1.36 15 289 454.68 12.29 0.65 296.07 -7.07 8 618 1.68 1.68 1.71 16 293 459.11 11.51 0.68 319.97 -26.97 9 770 2.03 2.01 2.05 17 279 474.43 11.89 0.57 268.13 10.87 10 564 1.45 1.21 1.24 18 552 483.64 11.62 1.16 567.61 -15.61 11 327 0.82 0.90 0.91 19 674 495.61 11.66 1.36 671.60 2.40 12 235 0.57 0.57 0.58 20 827 502.27 11.16 1.71 869.82 -42.82 13 189 0.45 0.98 21 1000 508.08 10.62 2.05 1054.96 -54.96 14 326 0.76 22 502 494.82 8.23 1.23 652.10 -150.10 15 289 0.65 23 512 515.48 9.48 0.91 455.74 56.26 16 293 0.65 24 300 523.76 9.36 0.58 303.43 -3.43 17 279 0.60 t Dự báo 25 359 577.75 13.82 0.49 253.08 105.92 18 552 1.17 49 396.25 26 264 559.18 10.58 0.60 363.51 -99.51 19 674 1.39 50 476.33 27 315 553.21 8.93 0.64 368.52 -53.52 20 827 1.67 51 525.10 28 361 555.81 8.29 0.68 382.52 -21.52 21 1000 1.98 52 578.44 29 414 596.14 11.50 0.58 322.43 91.57 22 502 0.97 30 647 597.22 10.46 1.16 707.64 -60.64 23 512 0.97 31 836 609.40 10.63 1.36 824.28 11.72 24 300 0.56 32 901 601.54 8.78 1.69 1058.95 -157.95 xt f 33 1104 596.16 7.36 2.03 1248.75 -144.75 1600 34 874 624.60 9.47 1.25 744.10 129.90 1400 35 683 656.60 11.72 0.93 579.97 103.03 1200 36 352 656.57 10.55 0.57 385.95 -33.95 1000 37 332 669.37 10.77 0.49 326.47 5.53 800 38 244 625.42 5.30 0.58 408.26 -164.26 600 39 320 604.72 2.70 0.63 403.07 -83.07 400 200 40 437 614.96 3.45 0.68 411.46 25.54 0 0 102030405060 41 544 681.07 9.72 0.61 361.08 182.92 42 830 696.18 10.26 1.16 798.88 31.12 Series1 300.00 43 1011 714.05 11.02 1.36 959.33 51.67 200.00 44 1081 708.22 9.34 1.67 1223.10 -142.10 100.00 45 1400 712.20 8.80 2.02 1454.23 -54.23 0.00 46 1123 756.55 12.35 1.27 900.92 222.08 0 102030405060 -1 0 0 .0 0 47 713 768.91 12.36 0.93 712.94 0.06 -2 0 0 .0 0 48 487 794.89 13.72 0.58 447.95 21
  25. Nhận xét: Sai số là các đại lượng ngẫu nhiên, có biên độ tăng dần. Nguyên nhân là do số quan sát dùng để tối ưu các tham số α, β, γ là quá ít (chỉ có 2 thời vụ) 2.5 Các phương pháp phân ly (Decomposition) 2.5.1 Các công thức chung Các mô hình làm trơn đã xét trước đây đều dựa trên ý tưởng là nếu chuỗi thời gian có một định dạng (mẫu) thì mẫu này có thể được tách khỏi tính ngẫu nhiên bằng cách làm trơn các giá trị quá khứ. Tác dụng của việc làm trơn là loại bỏ thành phần ngẫu nhiên trong chuỗi rồi sử dụng mẫu cho việc dự báo. Các phương pháp làm trơn đều chưa nhận dạng được từng thành phần riêng biệt của mẫu. Trên thực tế, mẫu có thể được tách (phân ly) thành hai hoặc nhiều nhân tố, đặc biệt là khi xuất hiện các kiểu thời vụ trong dữ liệu. Trong nhiều tình huống, sẽ là rất tốt nếu người dự báo biết được tỉ lệ nào của dữ liệu tại thời điểm đã biết phản ánh mức tăng / giảm chung và tỉ lệ nào của dữ liệu chỉ đơn giản thể hiện sự dao động của thời vụ. Các phương pháp phân ly là một trong các cách dự báo cổ điển nhất. Các phương pháp này thường cố gắng nhận dạng 3 thành phần tách biệt của chuỗi thời gian là xu thế, chu kì và thời vụ. Xu thế là tính xuyên suốt của chuỗi như tăng, giảm, ổn định. Chu kì là thời kì tăng trưởng hay suy thoái của nền kinh tế, của một ngành công nghiệp; giai đoạn ElNino hay LaNina của khí hậu. Thời vụ là các dao động của các quan sát theo một chiều dài thời gian cố định (mùa, năm, ) Dựa trên giả thiết dữ liệu được cấu thành từ một mẫu cùng với sai số (ngẫu nhên) Dữ liệu = mẫu + sai số = hàm(xu thế, chu kì, thời vụ) + sai số Mô hình chung của các phương pháp phân ly là xt = f(Tt, Ct, St, Et) (2.15) Nhận xét: để nhận diện được thành phần chu kì, ta cần có ít nhất 10 năm số liệu. Trong dự báo ngắn hạn, thành phần xu thế Tt thường bao gồm luôn thành phần chu kì Ct. Dạng hàm chính xác của quan hệ (2.15) phụ thuộc vào phương pháp phân ly cụ thể được sử dụng. Ta có các mô hình sau đây i> xt = Tt + St + Et mô hình cộng tính ii> xt = Tt St Et mô hình nhân tính iii> xt = Tt St + Et mô hình nhân tính với sai số cộng tính Các mô hình nhân tính thường xuất hiện nhiều trong lĩnh vực kinh tế. Đối với mỗi loại mô hình trên, phải vẽ đồ thị để kiểm tra xem yếu tố thời vụ là cộng tính hay nhân tính. 22
  26. 2.5.2 Phương pháp phân ly cổ điển (Classical Decomposition). Phương pháp này phân ly chuỗi thời gian ra các thành phần như thời vụ, xu thế Dự báo được coi như ngoại suy tuyến tính của chuỗi thời gian trong qua khứ. Các bước tiến hành như sau: ¾ Tính các trung bình trượt trung tâm Mt của chuỗi xt với độ dài N = L (độ dài của thời vụ) nhằm mục đích loại bỏ thành phần thời vụ và thành phần ngẫu nhiên, chỉ giữ lại thành phần xu thế (và chu kì). ¾ Nhận lại thành phần thời vụ và ngẫu nhiên: St + Et = xt – Mt đối với mô hình cộng tính thời vụ St + Et = xt / Mt đối với mô hình nhân tính thời vụ ¾ Cô lập các yếu tố thời vụ bằng cách lấy trung bình các yếu tố này tại các thời điểm cách nhau một khoảng thời gian L. Sau đó chuẩn hóa các yếu tố trung bình này để loại bỏ thành phần ngẫu nhiên có ở bước 2 (chỉ còn lại St, t=1, 2, , L). ¾ Tách thành phần thời vụ ra khỏi dữ liệu dt = xt - St nếu yếu tố thời vụ có dạng cộng tính dt = xt / St nếu yếu tố thời vụ có dạng nhân tính ¾ Tính thành phần xu thế Tt = a + bt bằng cách sử dụng phương pháp bình phương bé nhất tìm b và a từ hệ phương trình ⎧ N a+ (∑ t) b = ∑ d t ⎨ 2 ⎩(∑ t)a+ (∑∑ t )b = td t N t.d− t. d d t →b = ∑∑t ∑ t ; a =∑ t − b ∑ N∑∑ t2 − ( t)2 N N ¾ Tách thành phần xu thế Tt khỏi Mt để nhận được thành phần chu kì (nếu cần) Ct = Mt – Tt ¾ Tính sai số Et = Xt – Tt – St nếu yếu tố thời vụ có dạng cộng tính Et = Xt – Tt St nếu yếu tố thời vụ có dạng nhân tính ¾ Dự báo ft = Tt + St nếu yếu tố thời vụ có dạng cộng tính ft = Tt St nếu yếu tố thời vụ có dạng nhân tính 2.5.3 Các ví dụ Ví dụ 1: Sau khi phân tích chuỗi thời gian dưới đây, người ta nhận thấy chuỗi này có xu thế tuyến tính và yếu tố thời vụ cộng tính, L=3. Vậy mô hình phân ly thích hợp là xt = Tt + St + Et 23
  27. - Tính các trung bình trượt trung tâm Mt của xt-1, xt và xt+1, t = 2, 3, , 11 - Tính các yếu tố thời vụ St S1 = –2 là trung bình cộng của x1 – M1, x4 – M4, x7 – M7 và x10 – M10 S2 = 0 là trung bình cộng của x2 – M2, x5 – M5, x8 – M8 và x11 – M11 S3 = 2 là trung bình cộng của x3 – M3, x6 – M6, x9 – M9 và x12 – M12 → S1 + S2 + S3 = 0 có nghĩa là các St đã được chuẩn hóa sẵn! - Tách thành phần thời vụ ra khỏi dữ liệu bằng cách tính dt = xt - St - Để xác định thành phần xu thế Tt = a + bt ta sử dụng phương pháp bình phương bé nhất tính a và b: 2 ∑t = 78, ∑t = 650, ∑tdt = 1704.1, ∑dt = 217.5 → b = 2.02, a = 5.04 Từ đó ta lập được bảng tính dưới đây và tính được F13 = T13 + S13 = 5.04 +2.02 * 13 + (-2) = 29.26 F14 = T14 + S14 = 5.04 +2.02 * 14 – 0.03 = 33.24 t xt Mt xt-Mt St dt Tt Et 1 4.7 -2.00 6.7 7.06 -0.36 2 9.0 8.9 0.1 0.00 9.0 9.07 -0.07 3 13.0 11.1 1.9 2.00 11.0 11.09 -0.09 4 11.3 13.4 -2.1 -2.00 13.3 13.11 0.19 5 15.9 15.6 0.3 0.00 15.9 15.12 0.78 6 19.6 17.4 2.2 2.00 17.6 17.14 0.46 7 16.7 19.0 -2.3 -2.00 18.7 19.16 -0.46 8 20.7 20.8 -0.1 0.00 20.7 21.18 -0.48 9 25.0 23.1 1.9 2.00 23.0 23.19 -0.19 10 23.6 25.2 -1.6 -2.00 25.6 25.21 0.39 11 27.0 27.3 -0.3 0.00 27.0 27.23 -0.23 12 31.3 2.00 29.3 29.24 0.06 13 -2.00 31.26 Dự báo 29.26 a = 5.04 14 -0.03 33.28 33.24 b = 2.02 Chú ý: Khi tính trung bình trượt trung tâm Mt của chuỗi thời gian có độ dài thời vụ L là một số chẵn, xuất hiện vấn đề là đặt Mt vào đâu? Để khắc phục tình huống này, người ta lấy trung bình cộng của 2 trung bình trượt kề nhau làm giá trị của Mt. 24
  28. Ví dụ 2: Chuỗi thời gian dưới đây có yếu tố thời vụ với độ dài L=4. Hãy phân ly nó thành các thành phần, biết rằng nó có xu thế tuyến tính và thời vụ nhân tính. Ta áp dụng dạng mô hình phân ly: xt = Tt . St + Et với dt = xt / St Et = xt – Tt . St Lập bảng tính các thành phần, Mt là trung bình cộng của 2 trung bình trượt 4 bước kề nhau của các quan sát xt, t = 3, 4, , 14. 2 t xt 4-MA Mt xt/Mt St dt t tdt Tt et Dự báo 1 5 0.45 11.11 1 11.11 17.25 -2.76 7.76 2 21 21.25 1.19 17.59 4 35.18 19.09 -1.80 22.80 3 33 22.75 22 1.5 1.56 21.14 9 63.41 20.94 0.31 32.69 4 26 25.75 24.3 1.07 0.83 31.19 16 124.77 22.78 7.01 18.99 5 11 28.75 27.3 0.4 0.45 24.45 25 122.24 24.62 -0.08 11.08 6 33 27.25 28 1.18 1.19 27.64 36 165.84 26.46 1.40 31.60 7 45 28.25 27.8 1.62 1.56 28.82 49 201.75 28.31 0.80 44.20 8 20 31 29.6 0.68 0.83 24 64 191.96 30.15 -5.13 25.13 9 15 34.25 32.6 0.46 0.45 33.34 81 300.05 31.99 0.61 14.39 10 44 36.5 35.4 1.24 1.19 36.85 100 368.53 33.84 3.60 40.40 11 58 37.75 37.1 1.56 1.56 37.15 121 408.63 35.68 2.29 55.71 12 29 39.25 38.5 0.75 0.83 34.79 144 417.51 37.52 -2.28 31.28 13 20 43 41.1 0.49 0.45 44.45 169 577.88 39.37 2.29 17.71 14 50 43.25 43.1 1.16 1.19 41.88 196 586.30 41.21 0.80 49.20 15 73 1.56 46.76 225 701.34 43.05 5.78 67.22 16 30 0.83 35.99 256 575.88 44.89 -7.42 37.42 17 0.45 289 46.74 F17 = 21.03 136 497.1 1496 4852.39 18496 b= 1.84 a= 15.4 Nhận xét: Đồ thị của sai số et không có xu thế và không có định dạng (Et là các đại lượng ngẫu nhiên). 25
  29. 2.6 BÀI TẬP CHƯƠNG 2 1. Số máy tính bán ra hàng tuần của một đại lí cho bởi bảng dưới đây Tuần 1 2 3 4 5 6 7 8 9 10 11 12 Máy 75 74 79 83 69 78 71 80 77 85 81 70 Chứng tỏ rằng mô hình lảm trơn dạng mũ đơn là thích hợp cho dữ liệu này. Sử dụng mô hình làm trơn dạng mũ đơn với α = 0,1 dự báo số máy tính sẽ bán ra trong tuần tới. Có thể sử dụng mô hình này để cung cấp dự báo tin cậy cho tuần thứ 20 không? 2. Doanh thu bán hàng của một đại lí trong 6 tháng gần đây là 1000 3000 4500 5500 5100 5400 (tính bằng triệu đồng) Sử dụng đồ thị chứng tỏ rằng việc sử dụng tính xu thế là cần thiết Dự báo doanh thu cho 3 tháng tới, sử dụng lần lượt các phương pháp Trung bình trượt dạng mũ kép Mô hình Holt với α = 0,1 và β = 0,1 Mô hình Holt với α = 0,9 và β = 0,9 Tính các sai số MAE và MSE cho cả 3 phương pháp 3. Số lượng khách hàng đăng kí sử dụng truyền hình cáp trong 6 tháng đầu năm là 3403, 3790, 4025, 426, 4192, 4235. Hãy sử dụng mô hình Holt với các tham số trơn α = β =0,1 dự báo lượng khách hàng đăng kí sử dụ ng dịch vụ này trong tháng tới. Vẽ đồ thị minh họa các giá trị quan sát được và dự báo rồi rút ra nhận xét. Tính sai số MSE. Có cách nào cải thiện dự báo không? 4. Trung tâm báo chí quốc gia muốn dự báo nhu cầu hội nghị cho từng quý của năm tới. Dữ liệu thu thập được của 4 năm gần đây là 26
  30. Năm Quý Nhu cầuNăm Quý Nhu cầu 2002 1 10 2004 1 13 2 31 2 34 3 43 3 48 4 16 4 19 2003 1 11 2005 1 15 2 33 2 37 3 46 3 51 4 17 4 21 Vẽ đồ thị và giải thích vì sao mô hình Winters với xu thế tuyến tính, thời vụ cộng tính là thích hợp cho mục đích dự báo. Sử dụng mô hình này với các tham số trơn α = 0,2 ; β = 0,1 và γ = 0,1 dự báo nhu cầu hội nghị cho các quý của năm tới. 5. Sử dụng phương pháp phân ly cổ điển, tách chuỗi thời gian ở bài 4 ra các thành phần T, S, E. Vẽ đồ thị x ~ t để tìm mô hình phân ly thích hợp. Tính các dự báo cho 4 quý tiếp theo và đánh giá sai số. 6. Dữ liệu sau đây là tình hình tiêu thụ một loai sản phẩm theo quý của 6 năm gần đây Năm Quý Sản phẩm Năm Quý Sản phẩm Năm Quý Sản phẩm 1999 1 362 2001 1 473 2003 1 628 2 385 2 513 2 707 3 432 3 582 3 773 4 341 4 474 4 592 2000 1 382 2002 1 544 2004 1 627 2 409 2 582 2 725 3 498 3 681 3 854 4 387 4 557 4 661 Vẽ đồ thị để kiểm tra sự biến đổi theo thời vụ ? Biến đổi này là cộng tính hay nhân tính? Dữ liệu có thể hiện tính xu thế không ? Nếu có thì đó là xu thế tuyến tính hay dạng mũ? Sử dụng mô hình Winters thích hợp dự báo tình hình tiêu thụ sản phẩm cho 4 quý tới. 27
  31. 3 CHƯƠNG 3 : PHÂN TÍCH CHUỖI THỜI GIAN VÀ CÁC MÔ HÌNH CỦA BOX‐JENKINS Trong chương trước chúng ta đã xét các kỹ thuật dự báo dựa trên cơ sở các phép làm trơn với giả thiết rằng giá trị trung bình của chuỗi thời gian là hàm xác định của thời gian và quan sát ở bất kì thời điểm nào cũng là tổng của trị trung bình với thành phần sai số ngẫu nhiên xt = μt + et tức là sai số et là các biến ngẫu nhiên độc lập đối với t, còn μt là hàm trung bình xác định theo t. Nhận xét rằng nếu dãy {et} là các biến ngẫu nhiên độc lập thì dãy các quan sát {xt} cũng là các biến ngẫu nhiên độc lập. Nhưng giả thiết về sự độc lập của các quan sát {xt} lại thường không được đảm bảo vì có nhiều chuỗi thời gian mà các quan sát liên tiếp phụ thuộc vào nhau chặt chẽ. Trong các trường hợp đó, các kĩ thuật dự báo dựa trên các phép làm trơn có thể trở nên không thích hợp bởi chúng không tận dụng được ưu điểm của sự phụ thuộc giữa các quan sát một cách có hiệu quả nhất. Trên thực tế các phương pháp làm trơn thường chỉ cho kết quả tốt đối với các quan sát phụ thuộc vào thời vụ. Dựa trên giả thiết rằng các giá trị liên tiếp của chuỗi thời gian có liên quan với nhau, Box và Jenkins cố gắng khám phá điều đó và sử dụng trong các mô hình dự báo. Các mô hình này do Box và Jenkins đề xướng nên thường được gọi là các mô hình Box- Jenkins. Các kỹ thuật phân tích chuỗi thời gian xét trong chương này sẽ khai thác sự phụ thuộc giữa các quan sát. Các giá trị tương lai của chuỗi thời gian sẽ được xác định từ tổ hợp của các giá trị quá khứ và sai số quá khứ. 3.1 Các mô hình chuỗi thời gian ARMA (AutoRegressive‐Moving Average) Trong phân tích chuỗi thời gian, phương pháp Box – Jenkins, được đặt tên sau khi hai nhà thống kê học George Box và Gwilym Jenkins áp dụng các mô hình tự hồi quy trung bình trượt ARMA hay ARIMA, tìm ra mô hình phù hợp nhất của chuỗi thời gian các giá trị thời gian trong quá khứ để tạo ra các dự báo. 3.1.1 Mô hình tự hồi quy bậc p - AR(p) xt = Φ1 xt-1 + Φ2 xt-2 + +Φp xt-p + at + θ0 (3.1) Vậy mô hình AR(p) là tổng có trọng số của các quan sát đã cho với các trọng số (các tham số) Φ1, Φ2, , Φp. Các tham số này cần được ước lượng để tìm hàm dự báo. 28
  32. 3.1.2 Mô hình trung bình trượt bậc q - MA(q) xt = at – θ1 at-1 – θ2 at-2 – – θq at-q + θ0 (3.2) 3.1.3 Mô hình hỗn hợp tự hồi quy-trung bình trượt bậc (p,q) - ARMA(p,q) xt – Φ1 xt-1 – Φ2 xt-2 – – Φp xt-p = at – θ1 at-1 – θ2 at-2 – – θq at-q + θ0 (3.3) trong đó Xt là các giá trị quá khứ và hiện tại của chuỗi thời gian, θ0 là hằng số (mức của quá trình), Φt, θt là các hằng số (tham số) cần được ước lượng, 2 at là các nhiễu động (sai số) độc lập, có trung bình μA = 0, phương sai σA không đổi, và không nhất thiết phải có phân phối chuẩn. ¾ Một số mô hình thường gặp AR(1) xt = Φ1 xt-1 + at + θ0 AR(2) xt = Φ1 xt-1 + Φ2 xt-2 + at + θ0 MA(1) xt = at – θ1 at-1+ θ0 = (1– θ1B) at + θ0 MA(2) xt = at – θ1 at-1 – θ2 at-2 + θ0 2 =(1−θ1B −θ 2B )a t + θ 0 (mô hình trơn dạng mũ đơn) ARMA(0,0) xt = at + θ0 ARMA(1,1) xt = Φ1 xt-1 at - θ1 at-1+ θ0 ARMA(1,2) xt = Φ1 xt-1 + at - θ1 at-1 - θ2 at-2 + θ0 ARMA(2,1) xt = Φ1 xt-1 + Φ2 xt-2 + at - θ1 at-1 + θ0 3.2 Các điều kiện cần về tính dừng và tính khả nghịch Phương trình (3.3) có thể viết lại dưới dạng toán tử Φ(B) xt = θ(B) at + θ0 (3.4) 3.2.1 Điều kiện dừng 2 p Φ(B) = 1 – Φ1 B – Φ2 B – – Φp B được gọi là toán tử tự hồi quy, 2 B là toán tử lùi, tức là Bxt = xt-1; B xt = B(Bxt ) = Bxt-1 = xt-2, ; 29
  33. Box và Jenkins đã chứng minh rằng khi quá trình AR(p) là dừng (có trung bình và phương sai bất biến, các nghiệm của phương trình Φ(z) = 0 đều nằm ngoài vòng tròn đơn vị (đối với quá trình MA(q) điều kiện dừng luôn luôn thỏa mãn với mọi giá trị của các tham số θi). 3.2.2 Điều kiện khả nghịch 2 q θ(B) = 1 – θ1B – θ2B – – θqB được gọi là toán tử trung bình trượt; khi quá trình MA(q) là khả nghịch, các nghiệm của phương trình θ(z) = 0 đều nằm ngoài vòng tròn đơn vị (đối với quá trình AR(p), điều kiện khả nghịch luôn luôn thỏa mãn với mọi giá trị của các tham số Φi). Các điều kiện về tính dừng và tính khả nghịch được áp dụng một cách độc lập, và nói chung Φ(B) và θ(B) sẽ không có cùng bậc (p ≠ q). Chú ý: Với các điều kiện trên, các mô hình ARMA đều dừng theo nghĩa trung bình 2 2 μX= E(xt) và phương sai σX = V(xt) = E[(xt – E(xt)] là bất biến. Ví dụ: ● Xét mô hình ARMA(1,0) hay AR(1): xt = Φ1 xt-1 + at + θ0 Lấy kì vọng 2 vế Ext = Φ1 Ext-1 + Eat + Eθ0 hay μX – Φ1 μX = θ0 μX = θ0 /(1– Φ1) = const 2 Lấy phương sai 2 vế Vxt = Φ1 Vxt-1 + Vat + Vθ0 2 2 2 2 hay σX = Φ1 σX +σA + 0 2 2 2 σX = σA / (1 – Φ1 ) = const ● Xét mô hình MA(q): xt = at – θ1 at-1 – θ2 at-2 – – θq at-q + θ0 Lấy kì vọng 2 vế ta được μX = θ0 = const Lấy phương sai 2 vế ta được 2 = (1+ θ 2 + θ 2 + +θ 2) 2 = const σX 1 2 q σA ● Xét mô hình ARMA(1,1): xt = Φ1 xt-1 + at - θ1 at-1+ θ0 Lấy kì vọng 2 vế ta được μX = Φ1μX + θ0 → μX = θ0 /(1-Φ1) = const 2 2 2 2 2 2 Lấy phương sai 2 vế ta được σX = Φ1 σX +σA +θ1 σA 30
  34. (1+ θ2 ) σ2 → 2 1 A = const σX = 2 1− Φ1 → Từ nay về sau, để sử dụng được các mô hình ARMA, ta sẽ luôn luôn giả thiết chuỗi thời gian là dừng 3.3 Các trợ giúp cho việc phân tích chuỗi thời gian 3.3.1 Biểu diễn đồ họa chuỗi thời gian Vẽ đồ thị luôn luôn là khởi đầu cần thiết để xem xu thế, chu kì, thời vụ và các điểm ngoại lai của các quan sát. 3.3.2 Hệ số tự tương quan ACF (Auto Correlation Function) Định nghĩa 1: Gọi x1, x2, xN là tập các quan sát của chuỗi thời gian. Ta định nghĩa hệ số tự tương quan ρk = ρ( xt, xt + k) = E[xt - E(xt)][xt + k - E(xt + k)] / [σ(xt) σ(xt + k)] (3.5) = [E(xt xt + k) - E(xt) E (xt + k)] / [σ(xt) σ(xt + k)] với k = 1,2, (3.5’) Nhận xét: –1 ≤ ρk ≤ 1 với mọi k =1, 2, ρk ≈ 1: xt và xt + k có tương quan (dương) chặt; ρk ≈ –1: xt và xt + k có tương quan (âm) chặt; Đối với chuỗi thời gian dừng, trung bình μ và phương sai 2 là bất biến và do đó là X σX độc lập đối với t. 2 2 Vậy thì từ (3.5’) ta có ρ k = [E(xt xt + k) – μ X] / σ X với k = 1,2, (3.6) là các hệ số tự tương quan lí thuyết đối với các bước nhảy k = 1,2, Định nghĩa 2: Các hiệp phương sai (autocovariance) ở bước k là γk = Cov(xt, xt + k) = E[(xt – μX ) (xt + k – μX)] (3.6’) từ (3.5) ta có công thức nữa tính ρk ρk = γk / γ0 (cần nhấn mạnh rằng ρ0 = 1) (3.7) Một số ví dụ: k 2 2 ● Đối với quá trình AR(1) người ta tính được γk = Φ1 σA /(1– Φ1 ) k=0,1, k → ρk = Φ1 hay các hệ số tự tương quan của AR(1) tắt rất nhanh khi k >1 vì |Φ1| 1→ |Φ1| < 1 ). 2 2 ● Đối với quá trình MA(1) người ta tính được γ0 = σA (1+ θ1 ) 31
  35. ⎧ − θ ⎪ 1 , k= 1 → ρ = 2 ρ biến mất sau bước k = 1 k ⎨1+ θ1 → k ⎩⎪0, k> 1 ● Đối với quá trình ARMA(1,1) người ta tính được 2 γ0= Φ 1 γ 1 + σ A[1 −θ1( Φ 1 −θ 1)] 2 γ1= Φ 1 γ 0 − θ 1 σ A γk= Φ 1 γ k− 1 khi k≥ 2 Từ đó suy ra (1− Φ1 θ 1)( Φ 1− θ 1) ρ1 = 2 1+θ1 − 2 Φ1θ 1 2 k-1 ρk= Φ 1 ρ k− 1 = Φ 1 ρk− 2 = = Φ1 ρ1 khi k ≥ 2 Ước lượng của ρk là hàm tự tương quan mẫu SACF (Sample AutoCorrelation Function). Hệ số tự tương quan mẫu ở bước k được định nghĩa: N N− k 1 x ∑(xt − x)(xt+ k − x) ∑ t N− k t= 1 , k = 1,2, với⎯ x = t= 1 là trung bình mẫu rk = N 1 2 N ∑(xt − x) N r= 1 1 N− 1 1 N− 2 (xt − x)(xt+ 1 − x) ∑ (xt − x)(xt+ 2 − x) N− 1 ∑ N− 2 Cụ thể r1 = t= 1 , r2 = t= 1 , (3.8) N N 1 2 1 2 ∑ (xt − x) ∑ (xt − x) N r= 1 N r= 1 Ý nghĩa: rk đo mức độ tương quan giữa xt và xt+k. Vậy đối với chuỗi các biến ngẫu nhiên x1, x2, xN thì r1, r2, ≈ 0 k = 0 x1 x2 x3 x4 r0 = 1 x1 x2 x3 x4 k = 3 x1 x2 x3 x4 x5 r3 x1 x2 x3 x4 x5 Đối với chuỗi có yếu tố thời vụ theo quý: r3 là đáng kể 32
  36. Đối với chuỗi có yếu tố thời vụ theo năm: r12 là đáng kể Ví dụ: Chuỗi xt = {7 8 7 6 5 4 5 6 4} có các tự tương quan Một bước Hai bước Ba bước 7 8 7 7 7 6 8 7 8 6 8 5 7 6 7 5 7 4 6 5 6 4 6 5 5 4 5 5 5 6 4 5 4 6 4 4 5 6 5 4 6 4 r = 0.62 r = 0.32 r = 0.15 1 2 3 3.3.3 Hàm tự tương quan riêng phần PACF (Partial AutoCorrelation Function) Định nghĩa 3: Hàm tự tương quan riêng phần ρmm (bậc m) là hệ số tự tương quan cuối cùng Φm của mô hình AR(m). Ví dụ: AR(1) xt = Φ1 xt-1 + at → ρ11 = Φ1 AR(2) xt = Φ1 xt-1 + Φ2 xt-2 + at → ρ22 = Φ2 AR(m) xt = Φ1 xt-1 + Φ2 xt-2+ + Φm xt-m + at → ρmm = Φm • Có thể tính được Φ1, Φ2, , Φm, ví dụ đối với AR(1), Φ1 được tính như sau: - Nhân 2 vế của mô hình AR(1) với xt-1 ta được xt xt-1= Φ1 xt-1 xt-1 + at xt-1 - Lấy kì vọng 2 vế E(xt xt-1) = Φ1 E(xt-1 xt-1) + E(at xt-1) - Vì E(at) = 0 nên ta có phương trình tương đương γ1 = Φ1 γ0 → Φ1 = ρ11 = γ1 / γ0 = ρ1 → ta sử dụng xấp xỉ ρ11 = Φ1 ≈ r1 • Làm tương tự đối với mô hình AR(2) dẫn đến hệ phương trình Yule-Walker bậc 2 Φ1 + Φ2 ρ1 = ρ1 Φ1 ρ1+ Φ2 = ρ2 hay viết dưới dạng phương trình ma trận ⎛ 1 ρ1 ⎞ ⎛ Φ1 ⎞ ⎛ ρ1 ⎞ ⎜ ⎟ ⎜ ⎟ = ⎜ ⎟ ⎝ρ1 1 ⎠ ⎝ Φ 2 ⎠ ⎝ρ2 ⎠ 33
  37. giải ra ta được nghiệm Φ1 và Φ2 ρ()1− ρ Φ = 1 2 (bỏ qua) 1 2 1− ρ1 ρ− ρ2 Φ = 2 1 = ρ 2 2 22 1− ρ1 • Đối với mô hình AR(m): Giải hệ phương trình Yule-Walker bậc m ⎛ 1 ρ ρ ρ ⎞ ⎛ Φ ⎞ ⎛ ρ ⎞ ⎜ 1 2 L m− 1 ⎟ ⎜ 1 ⎟ ⎜ 1 ⎟ ⎜ ρ1 1 L ρm− 2 ⎟ ⎜ Φ 2 ⎟ ⎜ ρ2 ⎟ ⎜ ⎟ ⎜ ⎟ = ⎜ ⎟ và tính được Φm= ρmm ⎜ MMMMM ⎟ ⎜ MM⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ρm− 1 ρ m − 2 ρ m− 3 L 1 ⎠ ⎝Φ m ⎠ ⎝ρm ⎠ Ví dụ: Mô hình xt = 1.1 + 0.8 xt-1 + at Ý nghĩa: hệ số tự tương quan riêng phần ρmm dùng để đo độ tương quan giữa xt và xt+m khi các ảnh hưởng của m-1 bước giữa chúng bị loại bỏ. Trong việc phân tích chuỗi thời gian, ρmm giúp nhận dạng một mô hình ARMA thích hợp cho dự báo, chẳng hạn đối với mô hình AR(1) chỉ có Φ1 là đáng kể, AR(2) Φ1, Φ2 là đáng kể, MA(q) Φ1, Φ2, giảm về 0 theo dạng mũ. 34
  38. Nhận dạng mô hình AR(1) Nhận dạng mô hình AR(2) Trung Phương sai Tự tương quan Tự tương Mô hình 2 bình μX σX ρk quan riêng ρkk θ0 2 ρk = 0, k ≥ 1 ρkk = 0, k ≥ 1 ARMA(0,0) σA ρ11 = Φ1 2 2 k ARMA(1,0) θ0 /(1- Φ1) σA /(1- Φ1 ) ρk = Φ1 , k ≥ 1 ρkk = 0, k ≥ 1 2 ρ1 = Φ1/(1-Φ2) (1−Φ2) σ A ρ11 = Φ1/(1-Φ2) 2 2 (1+Φ2)[(1 −Φ 2) −Φ1] ARMA(2,0) θ0 /(1-Φ1-Φ2) ρk = Φ1 ρk-1+Φ2 ρk-2, k ≥ 1 ρ22 = Φ2 ρkk = 0, k>2 35
  39. 2 2 ⎧ − θ 2 (1+ θ1 ) σA 1 ρ11 =- θ1/(1+ θ 1 ) ρk = ⎪ 2 , k= 1 ⎨1 + θ1 ARMA(0,1) ⎪ θ0 ⎩0, k> 1 k 2 −θ1 (1 − θ1 ) ρkk = 2k+ 2 1− θ1 2 2 2 (1+θ1 +θ2 )σA − θ 1 (1 − θ 2 ) ρ11 = ρ1 ρ 1 = 2 2 1 +θ1 + θ 2 2 − θ 2 ρ2− ρ 1 ARMA(0,2) ρ 2 = 2 2 ρ = 1 +θ + θ 22 2 1 2 1 − ρ1 θ0 ρ k =0, k > 2 ρkk giảm dần về 0 (1 − Φ1 θ 1)( Φ 1− θ 1 ) ρ1 = 2 1 +θ1 − 2 Φ1θ 1 2 2 2 giảm dần về 0 ARMA(1,1) (1+θ1 -2Φ1θ 1) σ A ρk= Φ 1 ρ k− 1 = Φ 1 ρ k− 2 = 1−Φ2 = Φ k-1ρ khi k ≥ 2 θ0 /(1- Φ1) 1 1 1 2 Bảng các thông số μX, σX , ρk và ρkk của một số mô hình ARMA thường gặp: 3.3.4 Thống kê Q của Box-Pierce Giả sử mô hình ARMA(p, q) được sử dụng để mô phỏng chuỗi thời gian x1, x2, , Box và Pierce (1970) đã xây dựng một cách kiểm định giả thiết về sự phù hợp của mô hình: Các hệ số tự tương quan đầu tiên của phần dư đều bằng 0, tức là kiểm định giả thiết H0: ρ1 = ρ2 = = ρm = 0 m 2 bằng cách sử dụng thống kê Q= n∑ rk k= 1 với m = số bước lớn nhất được xét, n = N - d, N = số các quan sát ban đầu, d = số lần lấy sai phân chuỗi thời gian (sẽ xét sau), rk = hệ số tự tương quan mẫu ở bước k của phần dư {et= xt - ft }. 2 2 Q có phân phối χ với m-p-q bậc tự do (Q ~ χ m-p-q ) Nếu mô hình ARMA(p, q) được chọn là thích hợp với chuỗi thời gian {xt} thì giá trị Q tính được là nhỏ hơn giá trị tới hạn χα (α là mức ý nghĩa của kiểm định) và vì thế giả thiết H0 được chấp nhận, nếu trái lại, H0 sẽ bị loại bỏ Thống kê Q được sử dụng để kiểm định các hệ số tự tương quan của phần dư et sau khi tính được các dự báo ft. Nên sử dụng kiểm định này cùng với việc vẽ biểu đồ các hệ số rk của phần dư. 36
  40. 2 1 - α χ m – p – q chấp nhận giả thiết H0: đối thiết H1: không phải mọi hệ số ρkđều bằng 0 ρ1 = ρ2 = = ρm = 0 χα Sử dụng phân phối χ2 để kiểm định giả thiết về sự phù hợp của mô hình ARMA(p,q) 3.4 Các ứng dụng của các hệ số tự tương quan 3.4.1 Kiểm tra tính ngẫu nhiên của dữ liệu và phần dư ¾ Có thể chứng minh rằng 95% các hệ số tự tương quan mẫu rk nằm trong khoảng giới hạn (- 1.96 / n , 1.96 / n ) với n là số thời điểm quan sát. ¾ Khi một mô hình dự báo đã được chọn, các hệ số rk của phần dư (sai số) có thể tính được và do đó xét được tính ngẫu nhiên của chúng. Nếu chúng nằm trọn trong vùng giới hạn ±2/ n , tức là không đáng kể, thì các phần dư là ngẫu nhiên. Nếu các hệ số rk của phần dư có một định dạng xác định nào đó thì điều đó có nghĩa là mô hình dự báo đã chọn không phải là mô hình phù hợp. 3.4.2 Xác định tính dừng của chuỗi thời gian Các hệ số tự tương quan của chuỗi dừng sẽ giảm về 0 sau 2 hoặc 3 bước thời gian, trong khi đối với chuỗi không dừng, chúng vẫn còn khác 0 đáng kể sau một số bước. Trong việc phân tích chuỗi thời gian, nhận xét trên giúp người dự báo nhận ra một mô hình ARMA thích hợp. Chẳng hạn: MA(1) chỉ có ρ1 là đáng kể, MA(2) ρ1, ρ2 là đáng kể, AR(q) ρ1, ρ2, giảm về 0 theo dạng mũ, chỉ có ρkk, k = 1, , q là đáng kể ARMA(1, 1) cả ρk and ρkk đều tắt dần theo dạng mũ với k ≥ 2. 37
  41. Nhận dạng mô hình MA(1) Nhận dạng mô hình MA(3) 38
  42. Nhận dạng mô hình ARMA(1, 1) 3.4.3 Loại bỏ tính không dừng của chuỗi thời gian Bất cứ kiểu xu thế nào của các hệ số tự tương quan của {xt} đều do tính không xác thực của các tự tương quan chi phối. Nguyên nhân là do chuỗi thời gian đang xét chưa đảm bảo điều kiện dừng ban đầu (cả μ và 2 đều phải bất biến và do đó là độc X σX lập đối với t). Vì vậy phải loại bỏ tính không dừng của dữ liệu trước khi tiến hành phân tích bất kì chuỗi thời gian nào. Điều này được thực hiện thông qua việc sai phân hóa. Dấu hiệu nhận biết chuỗi thời gian không dừng là các hệ số rk của nó không tắt đi nhanh chóng sau vài bước đầu. Cần phải xác định bậc của sai phân đều và sai phân thời vụ để biến đổi chuỗi đã cho thành chuỗi mới có tính dừng. Trên thực tế, bậc của sai phân đều thường là d = 0, 1, hoặc 2 và chỉ cần xét đến 20 hệ số rk để kiểm tra kết quả. Chuỗi thời gian đã sai phân hóa thường được kí hiệu bởi Wt = Xt – Xt-1 khi d = 1 = Xt – B Xt = (1 – B) Xt (loại bỏ xu thế tuyến tính) Wt = (Xt+1 – Xt) – (Xt – Xt-1) khi d = 2 = Xt – 2Xt-1 + Xt-2 2 = (1 – B) Xt (loại bỏ xu thế bậc hai) k ở đây B là toán tử lùi (B Xt = Xt-1, B Xt = Xt-k) 39
  43. Nếu chuỗi thời gian không dừng do nguyên nhân phương sai 2 không bất biến, tức σX là tỉ lệ biến đổi của các quan sát không ổn định (được nhận ra khi vẽ đồ thị dữ liệu) thì cần thực hiện một phép biến đổi. Các phép biến đổi thường là yt = log xt α yt = x t , với α = 0,5 hoặc 0,25 Nếu có thể được, các phép biến đổi này nên thực hiện trước khi sai phân hóa. Các dự báo nhận được từ chuỗi đã biến đổi sẽ được chuyển lại (sử dụng phép biến đổi ngược) thành các dự báo cho chuỗi ban đầu. 3.4.4 Nhận biết tính thời vụ trong chuỗi thời gian Yếu tố thời vụ là dạng tự lặp lại sau một số bước cố định của chuỗi thời gian. Đối với chuỗi dừng, có thể nhận ra tính thời vụ bởi dấu hiệu khác 0 đáng kể của các hệ số rk sau một số bước thời gian L ≥ 2 cố định. Bất cứ hệ số rk nào khác không đáng kể đều chứng tỏ sự tồn tại của một định dạng dữ liệu. Để loại bỏ tính không dừng trong chuỗi thời gian chứa yếu tố thời vụ, ta cần sử dụng phép sai phân thời vụ. Chẳng hạn 12 D phép toán (1 - B ) sẽ loại bỏ các thành phần năm r12. Trên thực tế, bậc của sai phân thời vụ D = 1 và chỉ cần xét đến 36 hệ số rk đầu để kiểm tra kết quả. Tránh lạm dụng phép sai phân thời vụ (D ≥ 2) vì điều đó sẽ làm mất các quan sát thực. Ví dụ: Độ dẻo của một loại sản phẩm cao su là một đặc tính quan trọng. Bảng dưới đây cho biết độ dẻo của 100 sản phẩm sau cùng. Hãy xây dựng một mô hình chuỗi thời gian ARMA phù hợp cho quá trình này (đọc từ trên xuống, từ trái sang phải) 29.33 30.80 32.43 33.61 28.17 19.98 30.45 32.44 36.54 28.58 25.76 36.61 29.39 35.70 30.76 29.00 31.40 23.45 33.68 30.62 31.03 30.83 23.52 29.29 20.84 32.68 33.22 28.12 25.12 16.57 33.56 30.15 29.94 27.23 25.23 27.50 27.28 30.56 30.61 31.79 26.75 33.66 32.30 29.06 32.52 30.55 36.58 31.58 28.48 30.28 28.94 29.04 27.99 32.01 26.14 28.50 28.08 24.13 31.89 19.03 28.19 30.28 29.20 31.72 24.34 26.13 29.35 34.30 29.02 31.53 27.79 33.60 26.41 31.92 31.95 27.63 30.29 28.78 24.28 31.68 29.89 20.11 21.28 22.69 29.10 28.18 17.51 21.71 26.60 23.15 26.65 23.71 21.47 28.86 26.74 30.01 24.22 24.71 28.27 32.44 40
  44. → x= 28,57 Nhìn vào đồ thị xt ~ t ta thấy chuỗi này không có tính xu thế và yếu tố thời vụ. xt ~ t xt 40.00 35.00 30.00 25.00 20.00 15.00 0 20406080100120 Các hàm tự tương quan mẫu và tự tưong quan riêng phần mẫu cho chuỗi các quan sát trên được tính trong bảng dưới đây: Các hệ số tự tương quan mẫu rk tính bởi SIBYL Các hệ số tự tương quan riêng phần mẫu rkk tính bởi SIBYL Từ kết quả trên có thể đoán nhận mô hình cho chuỗi này là AR(2) và tính được r1 = 0,495 ≈ Φ1 và r2 = - 0,048 Giải hệ phương trình Yule-Walker để nhận được các ước lượng cho Φ1 và Φ2: 41
  45. Φ1 + Φ2 ρ1 = ρ1 Φ1 + 0,495 Φ2 = 0,495 Φ1 ≈ 0,6871 Φ1 ρ1+ Φ2 = ρ2 0.495 Φ1 + Φ2 = - 0,048 Φ2 ≈ –0,388 Ước lượng ban đầu của θ0 được suy ra từ công thức μX = θ0 / (1 – Φ1 – Φ2) → θ0 ≈ x (1 – Φ1 – Φ2) = 28,57(1 – 0,7 + 0,4) = 20 Ta có mô hình AR(2): xt = 20 + 0,495xt-1 – 0,388xt-2 + at Mô hình này có thể dùng để phát ra các giá trị ft khớp với chuỗi thời gian và tính phần dư at. Nhìn vào biểu đồ các hệ số tự tương quan mẫu rk của phần dư at tính bằng SIBYL ta thấy chúng đều không đáng kể và do đó at là các đại lượng ngẫu nhiên. Có thể xét thống kê Q của 10 hệ số tự tương quan đầu tiên: m 2 2 2 2 Q= n∑ rk = 100( 0 , 495+ 0 , 048 + +0 , 108 ) = 100 – 0,509714 = 50,97 k= 1 2 So sánh Q = 50,97 với giá trị tới hạn χ 0,05; 10-2 = 15,51 ta có thể kết luận mô hình xt = 20 + 0,495xt-1 – 0,388xt-2 + at là chưa phù hợp với chuỗi thời gian đã cho. Các hệ số tự tương quan mẫu rk của phần dư tính bởi SIBYL Nhận xét: Các tham số θ0, Φ1, Φ2 của mô hình AR(2) cho chuỗi thời gian này có thể cải thiện được bằng phương pháp bình phương bé nhất (tuyến tính). Để ý rằng xt = θ0 + Φ1xt-1 + Φ2xt-2 + at hay at = xt – θ0 – Φ1xt-1 – Φ2xt-2 là mô hình tuyến tính đối với các thông số θ0, Φ1, Φ2 → Tìm các ước lượng cho θ0, Φ1, Φ2 sao cho tổng các bình phương sai số 100 100 2 2 SSE=∑ at = ∑ (xt −θ 0 −ΦΦ 1x t− 1 − 2 x t− 2 ) đạt giá trị nhỏ nhất t= 3 t= 3 → Hệ 3 phương trình đại số tuyến tính xác định các tham số θ0, Φ1, Φ2 làm cực tiểu hàm SSE là 42
  46. ⎧ nθ0+Φ 1∑x t− 1 +Φ 2∑x t− 2 = ∑ x t ⎪ 2 ⎨θ0∑∑x t− 1+Φ 1x t− 1 +Φ2∑x t− 1 x t − 2= ∑ x t− 1 x t ⎪ 2 ⎩θ0∑∑x t2− +Φ 1x t1t2− x − +Φ 2∑x t2− = ∑ xt− 2 x t ⎧ 98θ0 + 2795, 20Φ1 + 2797, 79Φ2 = 2807, 66 ⎪ ⎨2795, 20θ0 + 8140116, Φ1 + 80643, 65Φ2 = 80925, 08 ⎪ ⎩2797, 79θ0 + 80643, 65Φ1 + 81546, 38Φ2 = 80074, 25 Giải hệ này ta tìm được θ0 = 20,64 Φ1 = 0,67 Φ2 = –0,38 là bộ thông số tối ưu 3.5 Các mô hình ARIMA 3.5.1 Các mô hình ARIMA không có tính thời vụ Dạng mô hình chung nhất mô phỏng quá trình hỗn hợp (cả AR và MA) là wt – Φ1 wt-1 – Φ2 wt-2 – – Φp wt-p = at – θ1 at-1 – θ2 at-2 – – θq at-q+ θ0 (3.7) hay (dạng toán tử) p q (1 – Φ1B – – ΦpB )wt = (1 – θ1B – – θqB ) at + θ0 (3.8) được gọi là mô hình ARIMA bậc (p,d,q) hay ARIMA(p, d, q), với d wt = (1–B) Xt p = bậc của quá trình AR d = bậc của sai phân đều để {wt} là chuỗi dừng q = bậc của quá trình MA I: tích hợp (Intergrated), vì {xt} đã được sai phân hóa để tạo ra chuỗi dừng. θ0 là hằng số, sẽ bằng 0 nếu và chỉ nếu trung bình của {wt} bằng 0. Nói cách khác: Các mô hình chứa xt (d = 0) gọi là các mô hình ARMA. Các mô hình chứa wt (d > 0) gọi là các mô hình ARIMA. Ví dụ : Nhu cầu hàng tuần về nhựa (tính bằng tấn) của một công ty sản xuất dây cáp điện của 100 tuần liên tiếp được cho trong bảng dưới đây. Hãy xây dựng một mô hình chuỗi thời gian ARIMA phù hợp cho quá trình này x ~ t (1-B)x ~ t 8000 600 7500 400 7000 6500 200 6000 0 5500 0 20 40 60 80 100 120 5000 -200 4500 -400 4000 0 20 40 60 80 100 120 -600 43
  47. 5000 5657 6132 7411 4965 6010 6111 7233 4496 6109 5948 6958 4491 6052 6056 6950 4566 6391 6342 6927 4585 6798 6626 6814 4724 6740 6591 6757 4951 6778 6302 6765 4917 7005 6132 6870 4888 7045 5837 6954 5087 7279 5572 6551 5082 7367 5744 6022 5039 6934 6005 5974 5054 6506 6239 6052 4940 6374 6523 6033 4913 6066 6652 6030 4871 6102 6585 5944 4901 6204 6622 5543 4864 6138 6754 5416 4750 5938 6712 5571 4856 5781 6675 5571 4959 5813 6882 5627 5004 5811 7011 5679 5415 5818 7140 5455 5550 5982 7197 5443 Nhìn vào đồ thị ta thấy chuỗi quan sát có khuynh hướng lệch ở vài chỗ nên không có được trị trung bình ổn định. Nếu vẽ đồ thị biểu diễn hàm tự tương quan ta thấy nó giảm về 0 rất chậm. Vậy chuỗi thời gian này không dừng. Lấy sai phân bậc nhất chuỗi xt ta được chuỗi wt = (1– B)xt. Đồ thị biểu diễn chuỗi wt cho thấy đây là chuỗi dừng, không có tính thời vụ và có trung bình w= 4,47. 44
  48. Có thể coi hằng số θ0 ≈ 4,47. Các hàm tự tương quan mẫu và tự tưong quan riêng phần mẫu cho chuỗi wt trên được tính trong bảng dưới đây: Các hệ số tự tương quan mẫu rk của chuỗi wt tính bởi SIBYL Các hệ số tự tương quan riêng phần mẫu rkk của chuỗi wt tính bởi SIBYL Nhìn vào các đồ thị trên ta thấy các hệ số tự tương quan mẫu rk là không đáng kể sau bước k = 1 còn các hệ số tự tương quan riêng phần mẫu rkk giảm về 0 sau 4 bước đầu. Điều đó chỉ ra rằng có thể đoán nhận mô hình MA(1) là phù hợp cho chuỗi wt hay mô hình ARIMA(0,1,1) là phù hợp cho chuỗi xt wt = 4,47 + at – θ1at-1 hay xt = xt-1 + 4,47 + at – θ1at-1 2 Tham số θ được tính từ công thức ρ1 = –θ/(1+θ ) thay ρ1 ≈ r1 = 0,41 ta tìm được θ1 = –1,92 và θ2 = –0,52 ta chọn θ1 = –0,52 thỏa mãn điều kiện |θ| <1 của mô hình ARMA(1,1) và nhận được mô hình thử nghiệm wt = 4,47 + at + 0,52at-1 hay xt = xt-1 + 4,47 + at + 0,52at-1 Các hệ số tự tương quan mẫu rk của phần dư at tính bởi SIBYL 45
  49. Nhìn vào biểu đồ các hệ số tự tương quan mẫu rk của phần dư at tính bởi SIBYL ta thấy có 3 giá trị đáng kể là r6, r11và r19 do đó at chưa hẳn là các đại lượng ngẫu nhiên. Có thể xét thêm thống kê Q của 10 hệ số tự tương quan đầu tiên: m 2 2 2 2 Q= n∑ rk = (100 − 1)(0,106+ 0,019 + + 0,024 ) =99 . 0,1402 = 13,88 k= 1 2 So sánh Q = 13,88 với giá trị tới hạn χ 0,05; 10-1 = 16,92 ta có thể kết luận mô hình xt = xt-1 + 4,47 + at + 0,52 at-1 là phù hợp với chuỗi thời gian đã cho. Để cải thiện giá trị của tham số θ1, người ta sử dụng phương pháp bình phương bé nhất phi tuyến, tìm θ1 tối ưu theo mục tiêu làm cực tiểu hàm SSE. Nhận xét rằng từ mô hình MA(1) wt = 0,47 + at – θ1at-1, t = 2, 3, đặt zt = wt – 0,47 ta có mô hình zt = at – θ1at-1 hay at = zt + θ1at-1 Sai số dự báo được cho bởi phần dư at = zt – ft. Giả sử aˆ t là xấp xỉ của at, t ≤ n, phụ thuộc vào sự lựa chọn thông số θ1 và thỏa mãn phép đệ quy aˆ t = zt + θ1 aˆ t-1 với các dự báo lùi zt = 0 và aˆ t = 0 với mọi t ≤ 1 Æ aˆ1 = 0 aˆ2 = z2 + θ1 aˆ1 = z2 aˆ3 = z3 + θ1 aˆ 2 = z3 + θ1z2 aˆn = zn + θ1 aˆn-1 = zn + θ1(zn-1 + θ1zn-2) = 2 n−2 = zn + θ1zn-1 + θ1 zn-2 + + θ1 z2 * Vậy aˆ t là hàm đa thức của θ1. Khi đó thông số tối ưu θ1 là điểm làm cực tiểu hàm n 2 SSE = ∑ at theo phương pháp bình phương bé nhất phi tuyến. Giá trị ước lượng ban t =2 đầu θ1 = –0,52 sẽ được dùng làm điểm xuất phát trong quá trình dò tìm. 3.5.2 Các mô hình ARIMA có tính thời vụ Dạng mô hình thời vụ chung nhất với chu kì L = 12 mô phỏng quá trình hỗn hợp (cả AR và MA) gọi là mô hình ARIMA(p,d,q) (P,D,Q)12 p 12 12P (1 – Φ1B – – ΦpB ) (1– Λ1B – – ΛPB ) wt 46
  50. q 12 Q = (1 – θ1B – – θqB ) (1 – γ1B – – γQB ) at + θ0 (3.9) với d 12 D wt = (1 – B) (1 – B ) xt p, d, q và θ0 là tương tự như trên P = bậc của quá trình thời vụ AR D = bậc sai phân của quá trình thời vụ cần cho tính dừng Q = bậc của quá trình thời vụ MA I: tích hợp (Intergrated), vì {xt}đã được sai phân hóa để tạo ra chuỗi dừng Trên thực tế, các mô hình chung nhất là ARIMA(1,d,0) (1,D,0)12 , ARIMA(1,d,0) (0,D,1)12 , ARIMA(0,d,1) (1,D,0)12 , ARIMA(0,d,1) (0,D,1)12 với d và D nhận các giá trị 0 hoặc 1. Ví dụ 1: Cho dãy 48 số liệu một loại nước giải khát đóng chai bán ra hàng tháng (tính theo kiện) của một hãng trong 4 năm liền. Hãy xây dựng một mô hình chuỗi thời gian ARIMA phù hợp cho quá trình này. Tháng 1994 1995 1996 1997 1 143 189 359 332 2 138 326 264 244 3 195 289 315 320 4 225 293 361 437 5 175 279 414 544 6 389 552 647 830 7 454 674 836 1011 8 618 827 901 1081 9 770 1000 1104 1400 10 564 502 874 1123 11 327 512 683 713 12 235 300 352 487 47
  51. Kiện Kiện~tháng 1600 1400 1200 1000 800 600 400 200 0 tháng 0 102030405060 Hàm tự tương quan mẫu của chuỗi xt Từ đồ thị x ~ t và hàm tự tương quan mẫu của xt ta thấy chuỗi thời gian đã cho không dừng và có tính thời vụ với độ dài L = 12 do đó phép biến đổi 12 wt = (1 – B) (1 – B ) xt được sử dụng để tạo chuỗi dừng. Hàm tự tương quan mẫu của chuỗi wt Hàm tự tương quan riêng phần mẫu của chuỗi wt 48
  52. Từ hai đồ thị trên ta thấy cả rk lẫn rkk đều giảm về 0 theo k nên có thể chọn mô hình thử nghiệm của wt là ARMA(1,1), tức là (1 – Φ1 B) wt = (1 – θ1 B) at ↔ wt = Φ1 wt-1 + at – θ1 at-1 hay xt = (1 + Φ1) xt-1 – Φ1 xt-2 – (1 + Φ1) xt-13 + Φ1 xt-14 + at – θ1 at-1 Sử dụng dữ liệu đã cho ta tính được Φ1 = 0,94 và θ1 = 0,14 Hàm tự tương quan mẫu của phần dư at Nhìn vào biểu đồ các hệ số tự tương quan mẫu rk của phần dư at tính bởi SIBYL ta thấy chỉ có r3 chạm vào đường giới hạn cho phép, do đó có thể xem at là các đại lượng ngẫu nhiên. Ta có thể kết luận mô hình ARIMA(1,1,1) (0,1,0)12 là phù hợp với chuỗi thời gian xt đã cho. Ví dụ 2: Quãng đường bay hàng tháng(đơn vị tính là triệu km) của một hãng hàng không trong 7 năm liền được cho trong bảng dưới đây. Hãy tìm mô hình ARIMA phù hợp với chuỗi thời gian này. Tháng 1994 1995 1996 1997 1998 1999 2000 1 7269 8350 8186 8334 8639 9491 10840 2 6775 7829 7444 7899 8772 8919 10436 3 7819 8829 8484 9994 10894 11607 13589 4 8371 9948 9864 10078 10455 8852 13402 5 9069 10638 10252 10801 11179 12537 13103 6 10268 11253 12282 12953 10588 14759 14933 7 11030 11424 11637 12222 10794 13667 14147 8 10882 11391 11577 12246 12770 13731 14057 9 10333 10665 12417 13281 13812 15110 16234 10 9109 9396 9637 10366 10857 12185 12389 11 7685 7775 8094 8730 9290 10645 11594 12 7602 7933 8280 9614 10925 12161 12772 49
  53. x ~ t 18000 16000 14000 12000 10000 8000 6000 4000 2000 0 0 20406080100 Từ đồ thị x ~ t ta thấy chuỗi thời gian đã cho không dừng cả về trung bình và phương sai nên ta sử dụng phép biến đổi yt = lnxt để làm mất tính không dừng về độ lệch Hàm tự tương quan mẫu của chuỗi yt Tuy nhiên đồ thị hàm tự tương quan mẫu cho thấy chuỗi yt chưa phải là chuỗi dừng về trung bình, ta dùng phép biến đổi sai phân đều zt = yt – yt-1 Hàm tự tương quan mẫu của chuỗi zt Hàm tự tương quan mẫu của chuỗi zt vẫn giảm về 0 rất chậm do tính thời vụ, thể hiện ở độ lớn của tham số r12, r24, Điều đó cho thấy phép sai phân thời vụ 12 wt = (1- B ) zt = zt – zt-12 = (yt – yt-1) – (yt-12 – yt-13) là cần thiết để biến zt thành chuỗi dừng wt 50
  54. Hàm tự tương quan mẫu của wt Hàm tự tương quan riêng phần mẫu của wt Từ hai đồ thị trên ta thấy mô hình thử nghiệm cho chuỗi wt nên là ARIMA(0,0,1) (0,0,1)12 hay 12 wt = (1– θ1B) (1– γ12 B ) at Nói cách khác, mô hình thử nghiệm cho chuỗi yt nên là ARIMA(0,1,1) (0,1,1)12 12 12 ↔ (1– B) (1– B ) yt = (1– θ1B) (1– γ12 B ) at hay yt – yt-1 – yt-12 + yt-13 = at – θ1at-1 – γ12at-12 + θ1γ12 at-13 Sử dụng phương pháp bình phương bé nhất (phi tuyến), người ta tính được các hệ số của mô hình là θ1= 0,67 và γ12 = 0,56, tức là 12 wt = (1– 0,67B) (1– 0,56 B ) at = at – 0,67 at-1 – 0,56 at-12 + 0,38 at-13 hay yt – yt-1 – yt-12 + yt-13 = at – 0,67 at-1 – 0,56 at-12 + 0,38 at-13 51
  55. Hàm tự tương quan mẫu của phần dư at Nhìn vào biểu đồ các hệ số tự tương quan mẫu rk của phần dư at tính bởi SIBYL ta thấy chỉ có r9 vượt ra khỏi đường giới hạn cho phép, do đó at chưa hẳn là các đại lượng ngẫu nhiên. Xét thêm thống kê Q của 10 hệ số tự tương quan đầu tiên 10 2 2 2 2 x Q= n∑ rk =( 72 − 2 ) ( 0 , 098+ 0 , 085 + +0 , 139 ) = 72 0,1736 = 12,50 k= 1 2 So sánh Q = 12,5 với giá trị tới hạn χ 0,05; 10-2 = 15,51 ta có thể kết luận mô hình ARIMA(0,1,1) (0,1,1)12 là phù hợp với chuỗi yt = ln xt đã cho. 52
  56. 3.6 BÀI TẬP CHƯƠNG 3 1. Chứng minh rằng đối với mô hình AR(1) ta có σ2 2 A σX = 2 1− Φ1 và do đó điều kiện dừng của mô hình là –1 {xt} có phải là chuỗi dừng không? ii> Cho x20 = 35, hãy mô phỏng 5 giá trị tiếp theo của quá trình này biết at ~ N(0, 4) 3. Xét mô hình AR(1): xt = Φ1 xt-1 + at. Chứng minh rằng mô hình này tương đương với mô hình 2 3 xt = at + Φ1 at-1 + Φ1 at-2 + Φ1 at-3 + 4. Cho 20 quan sát của chuỗi thời gian (đọc theo hàng) 4.2 5.8 6.9 7.62 5.57 3.34 2 1.7 2.02 2.71 3.63 5.18 7.11 8.26 7.95 6.78 5.07 5.04 6.02 7.61 i> Tìm ước lượng Φ1 cho mô hình AR(1) xt = Φ1 xt-1 + at của chuỗi thời gian trên bằng cách cực tiểu tổng bình phương các phần dư at (ft = Φ1 xt-1) ii> Vẽ đồ thị phần dư at ~ t. Từ đó có thể kết luận at là các đại lượng ngẫu nhiên không? Tại sao? 5. Xét mô hình MA(1): xt = at – θ1 at-1. i> Chứng minh rằng mô hình này tương đương với mô hình 2 3 xt = at – θ1 xt-1 – θ1 xt-2 – θ1 xt-3 – ii> Tìm khoảng giá trị chấp nhận được của θ1 để mô hình này khả nghịch. iii> Tìm ước lượng của θ1 nếu biết r1 6. Hiển thị các mô hình sau đây theo nghĩa xt được biểu diễn theo các số hạng trước đó của nó và của các nhiễu động at ARIMA(0,0,0) ARIMA(1,0,0) ARIMA(0,0,1) ARIMA(0,0,2) ARIMA(1,0,1) ARIMA(1,1,0) ARIMA(0,1,1) ARIMA(1,1,1) ARIMA(2,1,0) ARIMA(2,1,1) ARIMA(1,1,2) 7. Doanh thu hàng quý của một tổng công ty trong 5 năm liền là (đọc theo hàng) 25758 26067 26390 26710 27023 27338 27657 27967 28290 28598 28922 29519 30209 30611 31021 31235 31551 32071 32268 32394 i> Dữ liệu đã cho có tính xu thế và yếu tố thời vụ không? Tại sao? ii> Sử dụng phép biến đổi wt = (1– B) xt Chuỗi {wt}có phải là chuỗi dừng không? 53
  57. iii> Tìm mô hình ARIMA phù hợp với chuỗi {xt} 8. Số khách du lịch đến Thái Lan hàng tháng qua đường hàng không (tính theo ngàn người) của 3 năm gần đây là Tháng 2000 2001 2002 1 112 115 145 2 118 126 150 3 132 141 178 4 129 135 163 5 121 125 172 6 135 149 178 7 148 170 199 8 148 170 199 9 136 158 184 10 119 133 162 11 104 114 146 12 118 140 168 i> Dữ liệu đã cho có tính xu thế và yếu tố thời vụ không? Tại sao? ii> Sử dụng phép biến đổi 12 wt = (1– B)(1– B ) xt iii> Chuỗi {wt} có phải là chuỗi dừng không? iv> Tìm mô hình ARIMA phù hợp với chuỗi xt 9. Số người (ngàn) truy cập vào trang web của một công ty trong 5 năm là Tháng 1996 1997 1998 1999 2000 1 9.56 41.47 72.66 62.61 96.30 2 12.48 46.14 71.25 69.07 96.09 3 13.64 52.62 65.48 77.36 99.27 4 18.80 59.01 62.68 80.39 104.77 5 25.04 60.20 56.60 85.28 105.51 6 30.33 58.53 49.90 84.44 105.19 7 34.08 56.98 49.82 86.59 109.16 8 40.10 57.82 51.87 88.05 110.78 9 42.40 60.50 57.74 90.83 115.77 10 41.36 63.29 58.24 93.05 122.75 11 39.25 66.55 58.31 94.65 126.85 12 38.20 68.65 59.91 96.66 132.57 Tìm mô hình ARIMA phù hợp với chuỗi {xt} 54
  58. 4 CHƯƠNG 4: CÁC PHƯƠNG PHÁP DỰ BÁO CỦA BOX‐JENKINS Đây là quá trình xây dựng mô hình bằng cách chọn ra một mô hình từ lớp các mô hình ARIMA. Kỹ thuật của Box-Jenkins là tiến trình xây dựng mô hình chứ không chỉ đơn thuần là tiến trình tìm kiếm mô hình phù hợp bởi vì các mô hình được xác định trên cơ sở dữ liệu chứ không phải trên cơ sở giả thiết. 4.1 Các khâu chính trong phương pháp Box‐Jenkins Bước 1- Đoán nhận thăm dò: dữ liệu quá khứ được sử dụng để nhận dạng thử một mô hình ARIMA thích hợp. Bước 2- Ước lượng: dữ liệu quá khứ được sử dụng để ước lượng các tham số của mô hình thử nghiệm. Bước 3- Kiểm tra dự đoán: các đánh giá khác nhau được dùng để kiểm tra sự thích hợp của mô hình thử nghiệm, và nếu cần thiết, gợi ý một mô hình tốt hơn rồi sau đó mô hình này lại được xem như một mô hình thử nghiệm mới. Bước 4- Dự báo: Khi đã chọn được mô hình cuối cùng, nó được sử dụng để dự báo các giá trị tương lai của chuỗi thời gian. Đoán nhận mô hình Ước lượng các tham số của mô hình Kiể m tra dự đoán Mô hình có thích hợp khôn g? Không Có Sử dụng mô hình tạo ra các dự báo 55
  59. Trên thực tế, nhiều chuỗi thời gian có thể được biểu diễn bằng những mô hình đơn giản. Mô hình với số tham số ít nhất thường được ưa chuộng hơn. Thông thường đối với các mô hình ARMA(p,q) ta chỉ cần xét p ≤ 2 và, hoặc q ≤ 2. Có thể cải thiện việc biểu diễn mô hình bằng cách sử dụng một phép biến đổi dữ liệu gốc phù hợp. Dữ liệu đãợ đư c biến đổi, nếu chưa có tính dừng, sẽ được sai phân hóa cho đến khi đạt được tính dừng bởi vì ta bắt buộc phải làm việc với chuỗi thời gian dừng. Các hệ số tự tương quan của chuỗi thời gian dừng sẽ tắt rất nhanh và không có cấu trúc nào cả tức là không có bất kì một dấu hiệu định dạng nào. 4.2 Các nguyên tắc lựa chọn mô hình ARIMA(p,d,q) phù hợp Phương pháp Box-Jenkins được xem là một trong những kỹ thuật có hiệu quả cao trong việc phát ra các dự báo chính xác và tin cậy. Sức mạnh của nó là ở chỗ nó đưa ra những thông tin giúp nhà phân tích chuỗi thời gian lựa chọn mô hình phù hợp với dữ liệu quan sát được. Đối với các phương pháp khác, nhà phân tích giả thiết một mô hình nào đó rồi tiến hành ước lượng các tham số của mô hình. Trong giai đoạn đầu tiên, ta nhận dạng một mô hình thử nghiệm bằng cách so sánh các hàm tự tương quan mẫu và tự tương quan riêng phần mẫu của chuỗi thời gian dừng với các hàm tự tương quan và tự tương quan riêng lí thuyết của các mô hình ARMA. Trong khi rút ra nhận xét về các hàm tự tương quan và tự tương quan riêng lí thuyết của các mô hình ARMA khác nhau, cần nhớ rằng i> Chuỗi thời gian {wt} được xét là dừng theo nghĩa có trung bình không đổi (thường là 0) và phương sai bất biến. 2 iii> Phần dư at thường là đại lượng ngẫu nhiên có phân phối chuẩn N(0,σA ) iii> Các at là độc lập, tức là cov(at, at-k) =0 → E(at. at-k) = 0 iv> Các at không có tương quan với các quan sát trước đó, tức là E(at.wt-k) = 0 với k > 0 v> Các wt có tương quan với các quan sát wt-k trước đó (do định nghĩa của các mô hình ARIMA) vi> Tham số tự tương quan lí thuyết ở bước k được định nghĩa là γ E(w w ) k t t− k k = 0, 1, 2, ρk = = 2 γ0 E(wt ) Dưới đây là bảng đồ thị của các hàm tự tương quan và tự tương quan riêng phần d lí thuyết đối với các mô hình ARIMA của chuỗi thời gian xt : wt = (1 – B) xt không có thành phần (P,D,Q)L (xt không có tính thời vụ ) ρk ρkk 56
  60. ARIMA(0,d,1) w = a –θ a a t t 1 t-1 t-1 ARIMA(0,d,2) wt=a–θ1at-1–θ2at-2 at-2 ARIMA(1,d,0) wt=Φ1wt-1 +at ARIMA(2,d,0)w =Φ w + Φ w +a t 1 t-1 2 t-2 t ρk ρkk 57
  61. ARIMA(1,d,1)w =Φ w +a –θ a t 1 t-1 t 1 t-1 ARMA(0,d,0) wt = at Tất cả các ρk đều bằng 0 (các wt là độc lập Tất cả các ρkk đều bằng 0 (các wt cách với nhau) nhau k bước là độc lập với nhau) d 12 D Một số mô hình ARIMA có tính thời vụ của chuỗi xt : wt = (1 – B) (1 – B ) xt Tên mô hình Dạng mô hình ARIMA(0, d, 0) (0, D, 0)12 wt = at 12 ARIMA(0, d, 0) (0, D, 1)12 wt = (1 – γ12 B ) at 12 24 ARIMA(0, d, 0) (0, D, 2)12 wt = (1 – γ12 B – γ24 B ) at 12 ARIMA(0, d, 0) (1, D, 0)12 (1 – Λ12 B ) wt = at 12 24 ARIMA(0, d, 0) (2, D, 0)12 (1 – Λ12 B – Λ24 B ) wt = at 12 12 ARIMA(0, d, 0) (1, D, 1)12 (1 – Λ12 B ) wt = (1– γ12 B ) at 12 12 P 12 12 Q ARIMA(0, d, 0) (P, D, Q)12 (1–Λ12 B – –Λ12 P B ) wt = (1– γ12 B -– – γ12 Q B ) a 12 ARIMA(0, d, 1) (0, D, 1)12 wt = (1– θ1 B) (1 – γ12 B ) at 12 ARIMA(1, d, 0) (1, D, 0)12 (1 – Φ1 B) (1– Λ12 B ) wt = at 12 ARIMA(1, d, 0) (0, D, 1)12 (1 – Φ1 B) wt = (1 – γ12 B ) at 12 ARIMA(0, d, 1) (1, D, 0)12 (1 – Λ12 B ) wt = (1 – θ1 B ) at 12 ARIMA(1, d, 1) (1, D, 0)12 (1 – Φ1 B) (1–Λ12 B ) wt = (1– θ1 B ) at 12 ARIMA(1, d, 1) (0, D, 1)12 (1 – Φ1 B) wt = (1– θ1 B ) (1 – γ12 B ) at p 12 12 P ARIMA(p, d, q) (P, D, Q)12 (1– Φ1 B – – Φp B ) (1– Λ12 B – – Λ12 P B ) wt q 12 12 Q = (1– θ1 B – – θq B ) (1– γ12 B – – γ12 Q B ) at Trên thực tế, các mô hình sau đây thường được sử dụng: ARIMA(0, d, 1) (0, D, 0)12 ARIMA(1, d, 0) (0, D, 0)12 ARIMA(1, d, 0) (0, D, 0)12 ARIMA(0, d, 1) (0, D, 0)12 ARIMA(0, d, 1) (0, D, 1)12 ARIMA(1, d, 0) (1, D, 0)12 ARIMA(1, d, 0) (0, D, 1)12 ARIMA(0, d, 1) (1, D, 0)12 với d và D nhận các giá trị 0 hoặc 1. 4.3 Các hàm dự báo của các mô hình ARMA(p,q) Sau khi tìm được mô hình ARIMA phù hợp (tối ưu theo nghĩa làm cực tiểu hàm MSE), mô hình này sẽ được sử dụng để phát ra các dự báo cho các quan sát tương lai. 58
  62. Kí hiệu T là thời điểm hiện tại và T + τ là thời điểm cần dự báo trong tương lai. Gọi fT(τ) là dự báo cho thời điểm T+τ. Giá trị fT(τ) sẽ được xây dựng thông qua các dự báo liên tiếp cho các bước T+1, T+2, , T+τ–1. Các dự báo này được tính bằng cách lấy kì vọng (có điều kiện) của chuỗi thời gian xt theo mô hình ARIMA viết ở thời điểm xuất phát T cho các bước thời gian T + i. Trong quá trình này Các quan sát xT+j chưa biết sẽ được thay thế bởi fT(i) Các phần dư đã xuất hiện aT- i = xT- i – fT- i-1(1) Các phần dư aT+i chưa biết sẽ được thay thế bởi 0. Ở thời điểm đầu tiên, ta giả thiết rằng aT- i = 0 với mọi T – i ≤ 0 (các phần dư ở các thời điểm 0, –1, –2, đều bằng 0). 4.3.1 Một số mô hình ARMA thường gặp: - ARMA(0,0) xt = θ0 + at ta có fT(1) = E[xT+1] = E[θ0 + aT+1 ] = μX fT(2) = E[xT+2] = E[θ0 + aT+2 ] = μX fT(τ) = E[xT+τ] = E[θ0 + aT+τ ] = μX (4.1) - ARMA(1,0) xt = θ0 + Φ1 xt-1 + at ta có fT(1) = E[xT+1] = E[θ0 +Φ1xT + aT+1 ] = θ0 + Φ1 xT = μX (1 – Φ1) + Φ1 xT = μX + Φ1 (xT – μX) fT(2) = E[xT+2] = E[θ0 + Φ1xT+1 + aT+2 ] = E[θ0 + Φ1fT(1) + aT+2] 2 = θ0 +Φ1 [μX + Φ1 (xT – μX)] = μX + Φ1 (xT – μX) τ fT(τ) = μX + Φ1 (xT – μX) (4.2) - ARMA(2,0) xt = θ0 + Φ1 xt-1 + Φ2 xt-2 + at ta có fT(1) = E[xT+1] = E[θ0 + Φ1xT + Φ2 xT-1+ aT+1 ] = θ0 +Φ1 xT + Φ2 xT-1 = μX (1 – Φ1 – Φ2 ) + Φ1 xT + Φ2 xT-1 = μX + Φ1 (xT – μX) + Φ2 (xT-1 – μX) (4.3) fT(2) = E[xT+2] = E[θ0 + Φ1xT+1+ Φ2 xT + aT+2] = E[θ0 + Φ1 fT(1) + Φ2 xT + aT+2] = μX (1 – Φ1 – Φ2) + Φ1 fT(1) + Φ2 xT = μX + Φ1 [fT(1) – μX] + Φ2 (xT – μX) (4.3’) fT(τ) = μX + Φ1 [fT(τ –1) – μX] + Φ2 [(fT(τ -2) – μX) ] (4.3”) - ARMA(0,1) xt = θ0 + at – θ1 at-1 Tính các at = xt – θ0+ θ1at-1 với a0 = 0 (và từ đó suy ra các at-1) Tính fT(1) ≈ xT+1 = θ0 + aT+1 – θ1aT = θ0 – θ1aT (4.4) 59
  63. f100(τ) ≈ x100 + k = θ0 = μX ∀τ > 1(ở đây θ0 = μX) (4.4’) - ARMA(0,2) xt = θ0 + at – θ1 at-1 – θ2 at-2 ở đây θ0 = μX Tính các at với a0 = 0, at = xt – θ0 + θ1at-1 + θ1at-1 ta có fT(1) = E[xT+1] = E[θ0 + aT+1 – θ1 aT – θ2 aT-1] = μX – θ1 aT – θ2 aT-1 (4.5) fT(2) = E[xT+2] = E[θ0 + aT+2 – θ1 aT+1 – θ2 aT] = μX – θ2 aT (4.5’) fT(τ) = θ0 = μX với mọi τ > 2 (4.5”) - ARMA(1,1) xt = θ0 + Φ1 xt-1 + at – θ1 at-1 Tính các at với a0 = 0, at= xt – Φ1 xt-1 – θ0 + θ1at-1 + θ1at-1 ở đây θ0 = μX (1 – Φ1) ta có fT(1) = E[xT+1] = E[θ0 + Φ1 xT + aT+1 – θ1 aT ] = θ0 + Φ1 xT – θ1 aT = μX (1 – Φ1) + Φ1 xT – θ1 aT = μX + Φ1 (xT – μX) – θ1aT (4.6) fT(2) = E[xT+2] = E[θ0 + Φ1 xT+1 + aT+2 – θ1 aT+1 ] = θ0 + Φ1 xT+1 = μX (1 – Φ1) + Φ1 fT(1) = μX + Φ1 (fT(1) – μX) fT(τ) = μX + Φ1 (fT(τ-1) – μX) với mọi τ ≥ 2 (4.6’) 4.3.2 Giới hạn cho phép của các dự báo Một trong các cách tính các giới hạn cho phép của dự báo là viết lại mô hình đã cho dưới dạng tổng theo trọng số xT+τ = aT+τ + Ψ1 aT+τ-1 + Ψ2 aT+τ-2 + + Ψτ-1 aT+1 + Ψτ aT + Ψτ+1 aT-1 + (4.7) → fT(τ) = E[xT+τ] = Ψτ [xT – fT-1(1)] + Ψτ+1 [xT-1 – fT-2(1)] + (4.8) → eT(τ) = xT+τ – fT(τ) = aT+τ + Ψ1 aT+τ-1 + Ψ2 aT+τ-2 + + Ψτ-1 aT+1 4.4 Các ví dụ minh họa Ví dụ 1: Cho chuỗi các quan sát về độ dẻo của 100 sản phẩm có mô hình ARMA(2,0) phù hợp xt = 20,64 + 0,67 xt-1 – 0,38 xt-2 + at . Hãy dự báo độ dẻo cho hai sản phẩm tiếp theo. 20,64 Sử dụng công thức (4.3) và (4.3’) với T = 100, μX = = 28,71 1− 0,67 + 0,38 ta có f100(1) = μX + Φ1 (x100 – μX) + Φ2 (x99 – μX) = 28,71 + 0,67 (32,44 – 28,71 ) – 0,38 (26,74 – 28,71) = 31,95 f100(2) = μX + Φ1 [f100(1) – μX] + Φ2 (x100 – μX) 60
  64. = 28,71 + 0,67 (31,95 – 28,71 ) – 0,38 (32,44 – 28,71) = 29,44 Ví dụ 2: Cho các quan sát xt nhu cầu về nhựa hàng tuần (tính bằng tấn) của một công ty sản xuất dây cáp điện trong 100 tuần liên tiếp. Mô hình ARMA(0,1) phù hợp cho chuỗi wt = (1 – B) xt là wt = 4,47 + at + 0,52 at-1 (hay xt = xt-1 + 4,47 + at + 0,52 at-1) Hãy dự báo nhu cầu dùng nhựa cho hai tuần tới. Tính cột các at = wt + θ1at-1 với a0 = 0 và từ đó suy ra cột các at-1 Sử dụng các công thức (4.4) và (4.4’) ta có f100(1) = E[w100+1] = E[ θ0 + a101 – θ1 a100 ] = θ0 – θ1 a100 = 4,47 + 0,52 (–226,72) Æ x101 ≈ f100(1) + x100 = 5502,54 f100(2) = E[x100+2] = E[x100+1 + θ0 + a100+2 – θ1 a100 +1] = f100(1) + 4,47 Æ x102 ≈ θ0 + f100(1) = 5502,54+ 4,47 = 5507,01 Ví dụ 3: Cho 48 quan sát xt về số kiện một loại nước giải khát đóng chai bán ra hàng tháng của một hãng trong 4 năm liền. Hãy dự báo sự tiêu thụ loại nước này trong 3 tháng tới. 12 Mô hình ARMA(1,1) là phù hợp với chuỗi thời gian wt = (1–B)(1–B )xt. Nói cách khác, mô hình ARIMA(1,1,1) (0,1,0)12 là phù hợp với chuỗi thời gian xt. (1 – Φ1 B) wt = (1 – θ1 B) at ↔ wt = Φ1 wt-1 + at – θ1 at-1 (hay xt = (1 + Φ1) xt-1 – Φ1 xt-2 + xt-12 – (1 + Φ1) xt-13 + Φ1 xt-14 + at – θ1 at-1 ) với Φ1 = 0,94; θ0 = 0 và θ1 = 0,14 Tính cột các at = wt + θ1at–1 với a0 = 0 (và từ đó suy ra cột các at–1) Ở thời điểm dự báo T+τ ta có wT+τ = Φ1 wT+τ 1 + aT+τ – θ1 aT+τ 1 Với τ = 1 ta có fT(1) = E[wT+1] = Φ1wT – θ1 aT Vậy f48(1) = (0,94) w48 – 0,14 a48 = (0,94) (683) – (0,14 ) (–78,02) = 87,78 → x49 ≈ f48(1) + x48 + x37 – x36 = 554,78 Với τ = 2 ta có f48(2) = E[w48+2] = Φ1w48+1 = Φ1f48(1) = (0,94) (87,78) = 82,51 → x50 ≈ f48(2) + x49 + x38 – x37 = 549,29 61
  65. Ví dụ 4: Cho các quan sát xt về quãng đường bay hàng tháng (đơn vị tính là triệu km) của một hãng hàng không trong 7 năm liền. Mô hình ARIMA phù hợp với chuỗi thời gian yt = lnxt là 12 12 (1- B)(1- B ) yt = (1– θ1B)(1– γ12 B ) at với θ1 = –0,67 và γ12 = – 0,56 yt = yt-1 + yt-12 – yt-13 + at – 0,67 at-1 – 0,56 at-12 + 0,38 at-13 Hãy dự báo quãng đường bay của ba tháng tới. Ở thời điểm dự báo T + τ = 84 +1 ta có yT+1 = yT + yT-11– yT-12 + aT+1 – 0,67 aT – 0,56 aT-11 + 0,38 aT-12 → fT(1) = E[yT + yT-11 – yT-12 + aT+1 – 0,67 aT – 0,56 aT-11 + 0,38 aT-12] = yT + yT-11 – yT-12 – 0,67 aT – 0,56 aT-11 + 0,38 aT-12 Vậy f84(1) = y84 + y73 – y72 – 0,67 a84 – 0,56 a73 + 0,38 a72 = 9,46 + 9,29 – 9,41 – (0,67) (a84) – (0,56) (a73) + (0,38) ( a72) = 9,3485 Ở thời điểm dự báo T + τ = 84 +2 ta có yT+2 = yT+1 + yT-10– yT-11 + aT+2 – 0,67 aT+1 – 0,56 aT-10 + 0,38 aT-11 → fT(2) = E[yT+2] = fT(1) + yT-10 – yT-11 – 0,56 aT-10 + 0,38 aT-11 = 9,3485 + 9,2530 – 9,2910 –`(0,56) (0,0092) + (0,38) (0,0026) = 9,3063 Ở thời điểm dự báo T + τ = 84 +3 ta có yT+3 = yT+2 + yT-9 – yT-10 + aT+3 – 0,67 aT+2 – 0,56 aT-9 + 0,38 aT-10 → fT(3) = E[yT+3] = fT(2) + yT-9 – yT-10 – 0,56 aT-9 + 0,38 aT-10 = 9,3063 + 9,4716 – 9,2438 – (0,56) (0,0455) + (0,38) (0,0092) = 9,5484 Trở lại biến cũ ta nhận được các dự báo cho ba tháng tiếp theo là F84(1) = 11482 F84(2) = 11008 62
  66. F84(3) = 13128 63
  67. 4.5 BÀI TẬP CHƯƠNG 4 1. Cho các mô hình chuỗi thời gian xt phù hợp a) xt = 25 + 0,34 xt-1 + at Biết x100 = 28, tính các dự báo cho các thời điểm 101, 102, 103 b) xt = 15 + 0,36 xt-1 – 0,32 xt-2 + at Biết x50 = 32 và x51 = 30, tính các dự báo cho các thời điểm 51, 52 c) xt = 20 + at + 0,45 at-1 – 0,35 at-2 Biết x100 = 620 và x99 = 624, tính các dự báo cho các thời điểm 101, 102 2. a) Sau khi tìm được mô hình phù hợp với dữ liệu cho ở câu 7 chương 3, hãy phát dự báo cho 3 tháng sắp tới. b) Sau khi tìm được mô hình phù hợp với dữ liệu cho ở câu 8 chương 3, hãy phát dự báo cho 3 tháng sắp tới. c) Sau khi tìm được mô hình phù hợp với dữ liệu cho ở câu 9 chương 3, hãy phát dự báo cho 3 tháng sắp tới. 3. Doanh thu hàng tháng (đơn vị tính là triệu đồng) của một cửa hàng trong 8 năm liền được cho trong bảng dưới đây Tháng 1994 1995 1996 1997 1998 1999 2000 2001 1 26.0 29.0 27.0 27.5 23.9 27.6 26.9 28.4 2 24.5 24.7 26.3 27.2 24.7 23.4 26.6 25.5 3 27.9 31.3 29.8 30.2 27.5 25.0 27.6 26.6 4 29.1 32.4 32.6 28.6 26.7 26.0 27.1 26.2 5 34.7 33.9 35.1 34.1 28.7 31.0 29.8 29.3 6 33.1 35.0 34.4 30.9 30.3 29.3 29.1 28.8 7 36.0 36.4 35.7 34.7 31.3 31.7 32.6 31.2 8 37.5 36.5 33.6 33.7 32.1 32.0 31.6 31.9 9 34.8 34.4 31.9 33.6 31.2 30.0 31.1 28.5 10 35.5 33.9 35.1 31.0 31.4 31.8 33.2 32.2 11 33.4 33.9 33.4 28.9 30.8 33.6 33.6 31.7 12 32.9 36.4 37.6 29.7 30.6 31.6 34.0 31.8 Tìm mô hình ARIMA phù hợp với chuỗi thời gian này. Phát dự báo doanh thu cho 3 tháng sắp tới. 4. Số lượng khách (đơn vị ngàn người) của một hãng hàng không trong 3 năm là 112 118 132 129 121 135 148 136 119 109 104 111 115 126 141 135 125 149 170 170 158 133 114 141 145 150 178 163 172 178 199 199 189 162 146 161 Tìm mô hình ARIMA phù hợp với chuỗi thời gian này. Phát dự báo doanh thu cho 3 tháng sắp tới. 64
  68. 5 PHỤ LỤC: GIỚI THIỆU PHẦN MỀM DỰ BÁO SIBYL 5.1 Môi trường làm việc của Sibyl SIBYL/RUNNER là chương trình dự báo của Lincoln System Coporation gọn nhẹ, tương tác thân thiện với người sử dụng. Với Sibyl, ta có thể nhập dữ liệu mới từ bàn phím hay đọc từ file và ghi lại dữ liệu ra file với kiểu định sẵn là .SIB Có thể vào menu File để mở file dữ liệu có sẵn hay tạo file mới, hoặc vào menu Edit để soạn thảo dữ liệu. Sau đó chọn menu Analysis để xem các phân tích dạng dữ liệu (vẽ đường quan hệ x ~ t, phân tích yếu tố mùa, tính tự hồi quy, rồi vào Setup để thiết đặt các tính chất của dữ liệu và các tham số cần cho dự báo Menu Forecast sẽ hiển thị những mô hình dự báo có thể sử dụng cho file dữ liệu bạn đang mở để lựa chọn. Sau khi một mô hình dự báo đã được lựa chọn, Sibyl tự động tính toán và bạn mở menu Results để lựa chọn những kết quả muốn xem. 65
  69. Nếu các tiêu chuẩn dự báo chưa đạt mức đề ra, bạn có thể chỉnh sửa các thiết lập dữ liệu trong menu Setup, chọn mô hình dự báo khác và xem lại kết quả. 5.2 Một số phương pháp dự báo trong Sibyl 5.2.1 Các phương pháp trung bình trượt Trung bình trượt là lấy trung bình với số điểm được xét cố định. Trong chuỗi thời gian điều đó nghĩa là thông tin cũ nhất bị loại ra khi thông tin mới được cập nhật. Do đó giá trị trung bình nhận được luôn luôn là trung bình của một số không đổi các quan sát gần hiện tại nhất. ¾ Trọng số trung bìmh trượt là như nhau (1/N) cho n quan sát gần nhất. Vậy người sử dụng bắt buộc phải cung cấp đủ N quan sát gần bhất để dự báo bước tiếp theo. MAVE (Simple Moving Average). Phương pháp này thích hợp hơn cho chuỗi dữ liệu ổn định (không có tính xu thế) ¾ MAVE2 (Linear Brown's Moving Average). Đây là phương pháp trung bình trượt kép. Nó thích hợp cho chuỗi thời gian có tính xu thế. 5.2.2 Các phương pháp hồi quy tìm đường cong phù hợp với chuỗi dữ liệu (Trend-Cycle Regression Curve-Fitting Methods) ¾ SCURVE (Life Cycle Analysis). Phương pháp này giả thiết chuỗi thời gian có quan hệ dạng: x = ea + b/t 66
  70. với x là chuỗi thời gian, t là các thời điểm (1, 2, 3, , N), a và b là các hệ số tính được bằng phương pháp hồi quy Việc tìm đường cong phù hợp là thích hợp cho việc phân tích vòng đời của bất kỳ dạng chuỗi thời gian nào. Phương pháp này thường giả thiết rằng chuỗi dữ liệu xuất phát từ giá trị nhỏ nhưng lại thay đổi theo cấp số nhân, có một mức bão hòa. ¾ SREG (Straight-Line Trend Extrapolation). Phương pháp này giả thiết chuỗi thời gian có quan hệ dạng: x = a + bt với x là chuỗi thời gian, t là các thời điểm (1, 2, 3, , N), a và b là các hệ số tính được bằng cách lấy hồi quy ¾ EXGROW (Exponential Growth Trend Extrapolation). Phương pháp này giả thiết có một hằng số tỉ lệ tăng trong chuỗi thời gian và vì thế chuỗi thời gian có quan hệ dạng x = ea + bt với x là chuỗi thời gian, t là các thời điểm (1, 2, 3, , N), a và b là các hệ số tính được bằng phương pháp hồi quy Phương pháp này được dùng khi cần phải dự báo dài hạn một chuỗi thời gian có tỉ lệ tăng không đổi. 5.2.3 Các phương pháp làm trơn dạng mũ Các phương pháp này làm trơn (lấy trung bình) các giá trị quá khứ theo kiểu mũ, tức là chúng cho các quan sát ở xa hiện tại các trọng số giảm dần dạng mũ. ¾ EXPO (Single Exponential Smoothing) . Trong phương pháp này, giá trị dự báo được là trung bình trọng số của của một số giá trị đứng trước nó. Các trọng số này giảm dần theo dạng mũ khi các giá trị càng xa hiện tại. Công thức dự báo: f(t+1) = a x(t) + (1-a) f(t) với f(t+1) là dự báo, a là hằng trơn, t là thời điểm dự báo 67
  71. x(t) là giá trị của chuỗi tại thời điểm t. Phương pháp này chỉ tốt cho chuỗi thời gian ổn định, không có xu thế tăng / giảm. ¾ EXPO2 (Brown's 1-Parameter Linear Exponential Smoothing). Giống như EXPO, phương pháp này làm trơn (lấy trung bình) các giá trị quá khứ theo kiểu mũ, tức là chúng cho các quan sát ở xa hiện tại các trọng số giảm dần dạng mũ. Tuy nhiên nó thích hợp cho cả trường hợp dữ liệu có tính xu thế tuyến tính (bậc nhất). ¾ EXPOQ (Quadratic Exponential Smoothing). Đây cũng là phương pháp làm trơn dạng mũ nhưng ưu điểm của nó là cho phép dự báo đối với dữ liệu có xu thế là đường bậc hai. ¾ EXPOTL (Adaptive Response Rate Smoothing Trigg and Leach's exponential smoothing) là một phương pháp làm trơn dạng mũ với tỉ lệ thích hợp. Điểm khác biệt là nó không đòi hỏi người sử dụng phải chỉ ra một giá trị α (hoặc máy tính phải dò tìm một giá trị α tối ưu) phụ thuộc vào kiểu dữ liệu và biên độ biến đổi của nó. Vậy đây là một phương pháp tự thích nghi, không đòi hỏi thêm thông tin từ phía người sử dụng. ¾ EXPOW (Linear and Seasonal Winters' 3–Parameter Exponential Smoothing) Phương pháp này tương tự các phương pháp làm trơn dạng mũ kể trên, điểm khác biệt là nó xử lí được cả chuỗi dữ liệu vừa có nhân tố mùa vừa có tính xu thế. ¾ EXPOH (Linear Holt's Exponential Smoothing). Phương pháp Holt sử dụng nguyên tắc làm trơn dạng mũ đơn cho các quan sát quá khứ cho chuỗi dữ liệu có tính xu thế tuyến tính bằng cách sử dụng 2 tham số riêng biệt. ¾ EXPOD (Dampened-Trend Exponential Smoothing). Tương tự các mô hình của Winters và Holt là xem xét cả tính xu thế, nhưng nó thêm một tham số (là tổng của 3 tham số) để giảm thiểu tác động của xu thế khi dự báo dài hạn. Phương pháp này dùng cho trường hợp xu thế không có dạng đường thẳng mà là dạng dao động tắt dần theo thời gian. 5.2.4 Các phương pháp phân ly ¾ CENSUS (Census Decomposition). Phương pháp này tách chuỗi thời gian thành các thành phần như thời vụ, xu thế và ngẫu nhiên, nó cung cấp cả ước lượng cho các chỉ số về thời vụ. ¾ CENSUS II. Phương pháp này cung cấp các dữ liệu thống kê giúp điều chỉnh mô hình sau khi xem xét kết quả dự báo. Dự báo được coi như ngoại suy tuyến tính của chuỗi thời gian trong qua khứ ¾ DCOMP (Classical Decomposition). Phương pháp này phân ly chuỗi thời gian ra các thành phần như thời vụ, xu thế Dự báo được coi như ngoại suy tuyến tính của chuỗi thời gian trong qua khứ. 68
  72. 5.2.5 Phương pháp Box-Jenkins Box-Jenkins ARMA (Auto-Regressive Moving Average) là dạng chung nhất của các phương pháp dự báo chuỗi thời gian. Các giá trị tương lai của chuỗi thời gian được xác định từ tổ hợp của các giá trị quá khứ và sai số quá khứ. Dựa trên giả thiết rằng các giá trị liên tiếp của chuỗi thời gian có liên quan với nhau, Box-Jenkins cố gắng khám phá điều đó và sử dụng để dự báo. Phương pháp Box-Jenkins có thể dùng cho cả chuỗi thời gian có tính dừng cũng như không dừng. ¾ Cách sử dụng bộ công cụ ARIMA trong SIBYL Bước 1: Sử dụng menu File và Edit nhập dữ liệu và biển đổi nếu cần thiết (dựa trên minh họa đồ thị chuỗi quan sát và hàm tự tương quan mẫu) Bước 2: Mở menu Analyze để xem đồ thị cuỗi thời gian cũng như đồ thị hàm tự tương quan Bước 3: Mở menu Setup để thiết đặt các tính chất của dữ liệu như chọn chuỗi thời gian, số quan sát được sử dụng, xác định tính thời vụ, có phân li tính thời vụ hay không, số bước thời gian cần dự báo. 69
  73. Bước 4: Mở menu Forecast chọn phương pháp Box-Jenkins. Sử dụng chức năng AutoCorrelation để hiển thị các thông tin về tính dừng, tính thời vụ cho dạng mô hình ARIMA phù hợp, Chọn một mô hình thử nghiệm. 70
  74. Bước 5: Mở menu Results xem các thông số, các giá trị của mô hình khớp với các quan sát, các sai số, các giá trị của hàm tự tương quan mẫu và tự tương quan riêng phần mẫu và đồ thị minh họa để quyêt định lựa chọn mô hình này hay xét mô hình khác phù hợp hơn. 72
  75. Bước 6: So sánh các dự báo do các mô hình tạo ra để lựa chọn kết quả tốt nhất. 73