Luận văn Nghiên cứu ứng dụng mô hình arima để dự báo lượng mưa vụ đông xuân ở một số tỉnh vùng đồng bằng bắc bộ

h trượt bậc q (q = 1, 2,) và được ký hiệu là ARIMA(p,d,q). Dạng tổng quát của mô hình ARIMA(p,d,q) có thể được viết như sau [20]: Wt = µ + p1Wt-l + p2Wt-2 ++ ppWt-p - q1at-1 - q2at-2 -- qqat-q + at (2.1) Trong đó: Wt = Δ d yt d là bậc sai phân, μ là hằng số Với d = 0  Wt = yt ; với d = 1  Δyt = yt - yt-1; yt, yt-l, yt-2, , yt-p là giá trị quan trắc ở các bước thời gian t, t-1, t-2,, t-p at, at-1, at-2, , at-q là sai số ngẫu nghiên (giữa giá trị thực và giá trị tính toán) ở các bước thời gian t, t-1, t-2,, t-q; p1, p2, , pp ; q1, q3, ..., qq là các tham số hồi quy. Phương trình 2.1 cũng có thể viết gọn lại thông qua phép toán dịch chuyển lùi 23 tt a Bp Bq W )( )(   hoặc p(B)(1-B)d yt = μ + q(B)at (2.2) Trong đó: yt, at như đã trình bày ở trên B là phép tính dịch chuyển lùi: BWt = Wt-1 hay B k Wt = Wt-k p(B) = (1 – p1B – p2B 2 – – ppB p) là phép toán tự hồi quy q(B) = (1 – q1B – q2B 2 – – qqB q) là phép toán trung bình trượt 2.1.2. Mô hình động thái ARIMAX Mô hình động thái ARIMAX có sự khác biệt cơ bản so với mô hình tự hồi quy trung bình trượt ARIMA là ngoài việc xem xét quá trình tự hồi quy trung bình trượt của chuỗi yếu tố dự báo, nó còn cho phép xem xét ảnh hưởng của các chuỗi thời gian khác tác động đến yếu tố dự báo, chuỗi tác động (biến độc lập) được gọi là chuỗi nhập, chuỗi bị tác động (biến phụ thuộc) được gọi là chuỗi xuất. Giả sử ta có các chuỗi độc lập Xit (i = 1,2m; t=1,2n) và chuỗi phụ thuộc Yt (t=1,2n), khi đó mô hình động thái ARIMAX được viết dưới dạng tổng quát như sau: tti ki m i r i s i t a Bp Bq XB BS BU Y )( )( )( )( , 1     (2.3) Trong đó:  Yt là giá trị quan trắc ở các bước thời gian t; μ là hằng số;  B là phép toán dịch chuyển lùi theo quy tắc : BXt = Xt-1 , B k Xt = Xt-k ;  s isi1i0 B U B U U)( BU s i ; s isi1i0 BS BS S)( BS s i là những trọng số động thái của chuỗi độc lập thứ i;  k là thời điểm tác động của chuỗi độc lập thứ i tại thời điểm t = k;  p(B) = (1 – p1B – p2B 2 – – ppB p ); q(B) = (1 – q1B – q2B 2 – – qqB q ) là phép toán tự hồi quy và trung bình trượt của chuỗi phụ thuộc; 24  at, là sai số ngẫu nghiên (giữa giá trị thực và giá trị tính toán). Lưu ý : các chuỗi Xit và Yt trong công thức 2.3 phải là các chuỗi có tính ổn định ngẫu nhiên, nếu chuỗi không ổn định, sẽ cần phải thông qua bước sai phân để đưa chuỗi về dạng ổn định ngẫu nhiên. Tóm lại: Bản chất của các mô hình ARIMA và ARIMAX là mô hình ngẫu nhiên. Việc phân tích chuỗi thời gian trong các mô hình này bắt buộc phải chấp nhận một giả thiết hết sức cơ bản là tính ổn định của các quá trình ngẫu nhiên, tính ổn định ở đây có nghĩa là các đặc trưng thống kê (hay phân phối xác suất) không thay đổi theo thời gian. Trong thực tế nhiều quá trình ngẫu nhiên có tính ổn định trong một khoảng thời gian gián đoạn hữu hạn nào đó có thể coi là ổn định. Ví dụ chuỗi tổng lượng mưa tháng là chuỗi không dừng, còn chuỗi tổng lượng mưa năm có thể coi là dừng vì khi đó qui luật bên trong năm bị loại trừ. Các chuỗi không dừng có thể trở thành dừng nhờ một số phép biến đổi sai phân. Lợi thế cơ bản của các mô hình này là cho phép dự báo với độ chính xác nhất định, mặc dù chưa hiểu rõ bản chất của các quá trình tác động từ các nhân tố dự báo đến yếu tố dự báo. 2.2. Phƣơng pháp áp dụng mô hình ARIMA và ARIMAX đối với bài toán dự báo mƣa mùa Trong mục 2.1 đã trình bày các dạng tổng quát của mô hình ARIMA và ARIMAX, nó có thể bao gồm nhiều thành phần tham gia vào mô hình như: thành phần tự hồi quy, thành phần trung bình trượt, thành phần sai phân, thành phần ảnh hưởng của các chuỗi nhập khác (các chuỗi nhân tố dự báo), trong mỗi thành phần lại có các thành phần con khác nhau. Bài toán cần giải quyết ở đây là đưa ra được phương pháp xác định các thành phần có ý nghĩa về mặt thống kê để tham gia vào mô hình dự báo mưa hạn mùa. Đây là bài toán khá phức tạp, độ chính xác của mô hình dự báo không chỉ phụ thuộc vào các chuỗi nhập, chuỗi xuất mà còn phụ thuộc việc lựa chọn chính xác các thành phần tham gia vào mô hình dự báo. Để giải quyết bài toán này luận văn đã thực hiện theo các bước sau: 25 1) Áp dụng phương pháp thống kê sai phân để xác định tính ổn định ngẫu nhiên của các chuỗi dữ liệu tham gia vào mô hình ARIMA và ARIMAX; 2) Kế thừa phương pháp Box Jenkin đối với mô hình ARIMA và phương pháp Box Tao đối với mô hình ARIMAX trong việc nhận dạng các thành phần tự hồi quy, thành phần trung bình trượt và thành phần ảnh hưởng của các chuỗi nhập đến chuỗi lượng mưa thông qua việc xem xét sự biến đổi các hàm tự tương quan, tự tương quan riêng phần và tương quan chéo; 3) Sử dụng phương pháp bình phương tối thiểu trong việc xác định các tham số trong mô hình ARIMA và ARIMAX; 4) Áp dụng các phương pháp kiểm nghiệm giả thiết thống kê trong khí hậu để chọn lựa các tham số có đủ độ tin cậy thống kê tham gia trong mô hình ARIMA và ARIMAX; 5) Sử dụng công cụ phần mềm thống kê SAS để tính toán các đặc trưng của chuỗi thời gian và các tham số trong mô hình ARIMA và ARIMAX. Sau đây sẽ trình bày cụ thể từng nội dung này: 2.2.1. Xác định tính ổn định ngẫu nhiên của chuỗi thời gian Chuỗi thời gian là chuỗi số liệu được sắp xếp theo trình tự thời gian. Nếu một chuỗi thời gian có giá trị trung bình và phương sai không đổi theo thời gian thì chuỗi đó được xem là ổn định ngẫu nhiên (chuỗi có tính dừng) hay nói một cách khác cụ thể hơn đó là một chuỗi thời gian không có xu thế, không có chu kỳ, mà chỉ dao động xung quanh kỳ vọng của nó. Một chuỗi quan trắc khí hậu trung bình tháng thường bao gồm 3 thành phần: 1) thành phần ngẫu nhiên là sự tăng lên hay giảm đi thường xen kẽ nhau, góp phần làm cho các trị số khí hậu dao động xung quanh một giá trị nào đó. Giá trị đó có thể là trung bình số học, nếu chuỗi không có thành phần chu kỳ và xu thế. 2) Thành phần chu kỳ là những biến đổi của chuỗi lặp lại nhiều lần sau những khoảng thời gian nhất định nào đó. Mối tương quan giữa các thành phần trong một chu kỳ thường đạt trị số lớn nhất. 3) Thành phần xu thế là biểu hiện xu hướng tăng hoặc 26 giảm theo thời gian của các thành phần trong chuỗi, trị số đầu của xu thế là cực tiểu hoặc cực đại và trị số cuối của xu thế là cực đại hoặc cực tiểu. Biểu đồ minh họa 3 thành phần này được trình bày trong hình 2.1. Hình 2.1. Các thành phần trong chuỗi quan trắc khí hậu [10] Để loại bỏ thành phần xu thế và chu kỳ nhằm đưa các chuỗi quan trắc về dạng ổn định ngẫu nhiên, thường sử dụng phép lọc sai phân, phép lọc Loga, phép lọc căn thức...[10]. Trong luận văn này chúng tôi chọn phép lọc sai phân, cụ thể như sau: - Đối với việc loại bỏ thành phần xu thế: sử dụng phép biến đổi sai phân bậc 1 hoặc bậc 2. Sai phân bậc 1 là chênh lệch giữa 2 giá trị kề nhau trong chuỗi. ΔYt = Yt - Yt-1 (2.5) Trong đó: ΔYt là giá trị của sai phân bậc 1 Yt và Yt-1 là các thời đoạn trước và thời đoạn sau đó. a) b) c) d) 27 Nếu sai phân bậc 1 vẫn còn thể hiện xu thế thì thực hiện tiếp sai phân bậc 2. Sai phân bậc 2 chính là sai phân của sai phân bậc 1: Δ2(Yt) = ΔYt - ΔYt-1 = (Yt - Yt-1) - (Yt-1 - Yt-2) (2.6) Nếu sai phân bậc 2 chưa đạt được tính dừng ta có thể tiếp tục lấy sai phân bậc 3 hoặc cao hơn. - Đối với việc loại bỏ thành phần mùa và chu kỳ: Sai phân mùa là chênh lệch giá trị của hai quan trắc cách nhau khoảng thời gian L, L có thể là một năm, hai năm hay một mùa Ví dụ : nếu là số liệu tổng lượng mưa tháng, ta có L =12. Do đó sai phân mùa bậc 1 có tính mùa là: ΔYt = Yt - Yt-L = Yt - Yt-12 (2.7) Cũng có thể lấy sai phân bậc 2 của sai phân mùa bậc 1 khi chuỗi chưa đạt được độ ổn định: - Kiểm tra tính ổn định ngẫu nhiên của chuỗi Trong thực hành, để kiểm tra các chuỗi thời gian tham ra trong mô hình ARIMA hoặc ARIMAX đã đạt tiêu chuẩn ổn định ngẫu nhiên hay chưa, thường dựa vào hàm tự tương quan. Theo Quenouille đã chứng minh chuỗi thời gian được xem là ổn định ngẫu nhiên khi hầu hết hệ số tự tương quan của chuỗi (rk ) thỏa mản biểu thức giới hạn tin cậy (2.8) và tiến dần về 0, ngoại trừ một số bước trễ như bước mùa, vụ, chu kỳ , nằm ngoài khoảng này [6, 37]. Biểu thức giới hạn tin cậy có thể viết dưới dạng sau: rnkrn StrSt 1,2/1,2/    (2.8) n Sr 1  (2.9) Trong đó: Sr là sai số chuẩn của các hệ số tự tương quan rk ; 1,2/ nt là điểm phần trăm  =0.05 của phân bố Student với n-1 bậc tự do. Hình 2.2 minh họa chuỗi dữ liệu tổng lượng mưa tháng trước khi sai phân và sau khi sai phân. Trên hình này, phần hình A, A‟ ở phía trên ứng với trường hợp 28 chuỗi chưa ổn định, phần hình B, B‟ ứng với chuỗi sau khi sai phân và được xem là chuỗi ổn định ngẫu nhiên. Hình 2.2. Minh họa diễn biến của chuỗi lượng mưa tháng và hàm tự tương quan đối với trạm Hà Nội trước khí sai phân (A,A‟) và sau khi sai phân (B,B‟). 2.2.2. Nhận dạng cấu trúc của mô hình Sau khi đã loại bỏ được các thành phần chu kỳ, xu thế của chuỗi thời gian, sẽ tiến hành nhận dạng cấu trúc của mô hình. Box - Jenkin đã đưa ra phương pháp nhận dạng cấu trúc của mô hình ARIMA thông qua việc xem xét sự biến đổi của hàm tự tương quan (Autocorrelation function - ACF) và tự tương quan riêng phần (Part autocorrelation function - PAFC) để xác định các thành phần tự hồi quy (AR) và thành phần trung bình trượt (MA). Đối với mô hình ARIMAX, Box Tao đã đưa ra một số dáng điệu chính của hàm tương quan chéo (Cross correlation function - CCF) để xác định mức độ ảnh hưởng (hàm truyền) của các chuỗi nhập đến chuỗi yếu tố dự báo. Định nghĩa và thuật toán để tính các hàm ACF, PACF và CCF được trình bày chi tiết trong tài liệu ARIMA [37]. 29 Xuất phát từ bản chất của hàm ACF và PACF, Box - Jenkin đã đưa ra một số dạng biểu đồ thường gặp đối với hàm ACF và PACF, tương ứng với nó là các dạng của mô hình ARIMA nhằm hỗ trợ cho việc nhận dạng cấu trúc của mô hình, các dạng biểu đồ này được trình bày trong hình 2.3. Chi tiết về cách áp dụng biểu đồ này được trình bày dưới đây:  Nếu biểu đồ hàm ACF có dạng nhỏ dần theo các bước trễ thời gian và biểu đồ hàm PACF chỉ có giá trị khác 0 tại bước thời gian t-1, sau đó giảm đột ngột về 0, (cụm từ „khác 0‟ hay „bằng 0‟ ở đây được hiểu theo thuật ngữ thống kê, nếu các giá trị này nằm trong khoảng từ - 1,2/ nt *Sr đến + 1,2/ nt * Sr được xem là bằng 0, ngoài khoảng này được xem là khác 0 ) thì có một thông số tự hồi qui (p=1) được chọn, mô hình có dạng ARIMA(1,0,0). Ngược lại, khi biểu đồ hàm PACF tắt dần, hàm ACF có giá trị khác 0 bước t-1, sau đó giảm đột ngột về 0, trong trường hợp này mô hình có dạng ARIMA(0,0,1). Dáng điệu của hàm ACF và PACF trong các trường hợp này được minh họa trên hình (Hình 2.3a).  Tương tự như trên, nếu biểu đồ hàm ACF tắt dần, hàm PACF có giá trị khác 0 ở các bước thời gian t-1, t-2, sau đó giảm đột ngột về 0 thì mô hình có hai thông số tự hồi qui (p=2), mô hình có dạng ARIMA(2,0,0). Ngược lại, nếu hàm PACF có dạng tắt dần, hàm ACF có giá trị khác 0 ở các bước thời gian t-1, t-2, sau đó giảm đột ngột về 0 thì mô hình có hai thông số trung bình trượt (q=2), mô hình có dạng ARIMA(0,0,2). Đồ thị của hạn ACF và PACF trong các trường hợp này được minh họa trên hình (Hình 2.3b).  Khi biểu đồ hàm ACF có dạng tắt dần và có giá trị khác 0 ở các bước thời gian t-1, t-2t-p, tương tự hàm hàm PACF có dạng tắt dần và có giá trị khác 0 ở các bước thời gian t-1, t-2t-q, trong trường hợp này cả 2 thành phần AR và MA đều có trong mô hình, dạng của mô hình trong trường hợp này sẽ là ARIMA(p,0,q). Đồ thị của hạn ACF và PACF trong các trường hợp này được minh họa trên hình (Hình 2.3c). 30 (2.3a) (2.3b) (2.3c) Hình 2.3 Một số dạng chính của hàm ACF và PACF tưng ứng với các dạng mô hình ARIMA khác nhau [20] Đối với việc xác định ảnh hưởng của từng chuỗi nhập đến yếu tố dự báo trong mô hình ARIMAX, Box Tao cũng sử dụng phương pháp trực quan để xem xét sự biến đổi của hàm tương quan chéo (Cross correlation function - CCF), từ đó đưa 31 ra hàm truyền tương ứng của chuỗi nhập X tham gia trong mô hình ARIMAX. Nội dung của phương pháp có thể bao gồm 4 dạng chính sau: Hình 2.4. Một số dạng chính của hàm tương quan chéo giữa biến nhập (X) và biến phụ thuộc (Y) tưng ứng với các dạng mô hình ARIMA khác nhau [20] 1) Nếu hàm tương quan chéo (CCF) giữa biến độc lập X và biến phụ thuộc Y có giá trị „khác 0‟ tại bước thời gian (t-b), sau đó giảm đột ngột, các bước (2.4A) (2.4B) (2.4C) (2.4D) 32 thời gian khác đều có giá trị „bằng 0‟. Khi đó hàm truyền của biến X tham gia vào mô hình ARIMAX sẽ có dạng U0Xt-b, (hình 2.4A). 2) Nếu hàm CCF có giá trị „khác 0‟ tại bước thời gian (t-b) và (t-b-1), sau đó giảm đột ngột, các bước thời gian khác đều có giá trị „bằng 0‟. Khi đó hàm truyền của biến X tham gia vào mô hình ARIMAX sẽ có dạng (U0 + U1B)Xt-b, (hình 2.4B). 3) Nếu hàm CCF có giá trị „khác 0‟ tại bước thời gian (t-b), sau đó có xu thế giảm dần nhưng vẫn „khác 0‟ ở bước (t-b-1), còn các bước thời gian khác có giá trị „bằng 0‟. Khi đó hàm truyền của biến X tham gia vào mô hình ARIMAX sẽ có dạng sau : (hình 2.4C). 4) Nếu hàm CCF có giá trị „khác 0‟ và tại bước thời gian (t-b-1) và (t-b), đạt cao nhất tại (t-b-1), sau đó có xu thế giảm dần nhưng vẫn „khác 0‟ ở bước (t-b-2), còn các bước thời gian khác có giá trị „bằng 0‟. Khi đó hàm truyền của biến X tham gia vào mô hình ARIMAX sẽ có dạng sau : (hình 2.4D). 2.2.3. Xác định các tham số của mô hình Như đã trình bày ở trên (mục 2.2.2), tùy thuộc vào đặc tính và mối quan hệ của các chuỗi nhập và chuỗi xuất, mô hình ARIMA và mô hình ARIMAX có thể có một trong các dạng chính sau: Khi chỉ có thành phần AR(p), mô hình sẽ có dạng: yt = p1yt-l +p2yt-2 + +ppyt-p + at (2.10) Khi chỉ có thành phần MA(q), mô hình sẽ có dạng: yt = µ + q1at-1 + q2at-2 ++ qqat-q + at (2.11) Khi có đủ cả 2 thành phần AR và MA, mô hình sẽ có dạng: yt = µ - p1yt-l - p2yt-2 -- ppyt-p + q1at-1 + q2at-2 ++ qqat-q + at (2.12) btX BS BUU    1 10 1 btX BS U   1 0 1 33 Khi có đủ tất cả các thành phần AR, MA và thành phần động thái X1, X2,, Xm mô hình ARIMAX có thể có dạng: yt = µ - p1yt-l - p2yt-2 -- ppyt-p + q1at-1 + q2at-2 ++ qqat-q + (Ui1B+ Ui2B 2 ++ UihB h )Xit-b + at; (2.13) trong đó: (i = 1,2,,m; h = 1,2,3,). Việc xác định các tham số hồi quy trong các mô hình này được dựa theo nguyên tắc bình phương tối thiểu, trong đó các hệ số hồi quy được xác định sao cho tổng bình phương độ lệch giữa giá trị thực và giá trị mô phỏng là nhỏ nhất. Chi tiết về phương pháp này được trình bày trong sách giáo trình [10]. Tuy về nguyên tắc việc xác định các tham số trong mỗi thành phần của từng mô hình là giống nhau, nhưng về cách tính cụ thể có nhiều cách xử lý riêng. Do vậy ở đây sẽ trình bày tóm tắt cho một số dạng chính sau. 1) Xác định các thông số pi khi chỉ có thành phần AR(p) tham ra vào mô hình Đây là dạng hàm tương tự như hồi quy tuyến tính nhiều biến, Tuy nhiên các hệ số pi phải thỏa mản tiêu chuẩn hội tụ |p1 + pi ++ pi|<1 và được xác giải thông qua hệ phương trình Yule_Walker [6] Ck = p1Ck-l + p2Ck-2 + + ppCk-p; với k =1, 2,,p (2.14) Ck là các mô men tương quan (Covarian - hiệp phương sai) giữa biến phụ thuộc (yt) và các biến độc lập (yk-p). Ck-l , Ck-2 , Ck-p là các mô men tương quan giữa các biến độc lập với nhau. Đây là hệ p phương trình với các ẩn số là p1, p2 ,,pp và có thể giải được theo nhiều phương pháp khác nhau như: dùng phương pháp ma trận nghịch đảo, phương pháp khử Gauss, phương pháp Cramer 2) Xác định các thông số qi khi chỉ có thành phần MA(q) tham ra vào mô hình Các thông số qi thoả mãn một hệ phương trình tương tự như hệ Yule_Walker, được suy ra từ quan hệ: 22 2 2 1 11 0 k ...1 ... q qkpkk k qqq qqqqq r       (2.15) Trong đó : rk là hệ số tự tương quan của chuỗi yt, 34 γk và γ0 là các mô men tương quan bậc k và bậc 0 (k=0) của chuỗi Cho k = 1,2,...,q ta được một hệ phương trình phi tuyến. Phương pháp để giải hệ phương trình này được trình bày trong sách giáo trình [6,10]. 3) Xác định các thông số pi, qi khi cả 2 thành phần AR(p) và MA(q) tham ra vào mô hình - Đối với thành phần AR(p): để giải quyết độc lập các giá trị pi theo công thức truy hồi Durbin hay hệ phương trình Yule_Walker thì hệ thức (2.15) được viết bắt đầu từ k > q tức là khi ấy các bi = 0. Do đó ta có hệ : Cp+1 = p1Cp + p2Cp-2 + + ppCp+1-p Cp+2 = p1Cp+1 + p2Cp+ + + ppCp+2-p Cp+p = p1Cp+p-1 + p2Cp+p-2 + + ppCp (2.16) Hệ phương trình (2.26) cũng là hệ tuyến tính bậc nhất, do vậy có thể giải ra tìm các nghiệm p1, p2,, pp. - Đối với thành phần MA(q): để tìm các hệ số qi ta cũng xuất phát từ quan hệ (2.15), nhưng khác với mô hình MA(q), ở đây γk và γ0 (hay Ck và C0) không phải là của yt mà là của thành phần ngẫu nhiên at. Nghĩa là γk=γka, Ck=Cka. Quan hệ giữa γka và γky hay Cka và Cky có dạng [6]:     ip h kih p i ky p i ika daaCaC 010 2 (2.17) Với dk = Ck+i +Ck-i Cky chính là mô men tương quan bậc k của y và có thể tính theo công thưc: Cky = [(yt – ytb)(yt+k – ytb)] ; t = 1,2,...,n ; k = 1,2,,m ; (m < n) Sau khi có Cky, sẽ tính được Cka theo (2.17), thay vào (2.15) được một hệ phương trình, giải hệ này ta được các hệ số qi. 4) Xác định các thông số động thái Ui,t-b đối với các biến nhập Xi,t-b và các thông số pi, qi trong mô hình động thái ARIMAX 35 Đối với các thông số Ui,t-b tương ứng với các biến nhập Xi,t-b được xác định tương tự như mô hình hồi quy tuyến tính nhiều biến. Trong đó, biến phụ thuộc là chuỗi Yt, biến độc lập là các chuỗi Xi,t-b (i = 1,2m là ký hiệu biến nhập, t-b là độ chễ thời gian so với biến phụ thuộc). Phương pháp xác định các tham số hồi quy tuyến tính nhiều biến (Ui,t-b) được trình bày trong giáo trình [10]. Đối với việc xác định các thông số pi, qi trong mô hình động thái ARIMAX được xác định tương tự như trong mô hình ARIMA, tuy nhiên chuỗi tham gia vào mô hình ARIMA trong trường hợp này là phần sai số của mô hình hồi quy tuyến tính giữa Yt, biến độc lập Xi,t-b. 2.2.4. Kiểm định mô hình Việc kiểm định để đánh giá sự phù hợp của mô hình là vô cùng quan trọng, một mô hình được xem là phù hợp khi và chỉ khi thỏa mản các tiêu chuẩn như sau: 1) các giá trị sai số của mô hình phải độc lập nhau; 2) giá trị của mỗi thông số cần phải đủ lớn để mô hình có ý nghĩa thống kê; 3) các thông số của mô hình phải độc lập nhau. Dưới đây sẽ đi vào từng vấn đề cụ thể: 1) Kiểm tra tính độc lập của chuỗi sai số Việc kiểm tra tính độc lập của chuỗi sai số chính là xem xét hàm tự tương quan của các sai số đó. Nếu sai số có tính độc lập thì giữa chúng không có tương quan, hay nói cách khác các hệ số tương quan giữa chúng phải thỏa mản biểu thức giới hạn tin cậy (công thức 2.8 mục 2.2.1). Nếu không đạt tiêu chuẩn này thì giữa chúng có mối quan hệ với nhau, như vậy mô hình được chọn chưa phù hợp với các chuỗi số liệu được xem xét. Khi đó cần dựa vào hàm ACF và PACF của chuỗi sai số để điều chỉnh lại dạng của mô hình. 2) Kiểm định độ lớn của các thông số Thực chất các mô hình trong luận văn này là sử dụng phương pháp thống kê, vì vậy các thông số được xác định trong mô hình đều phải kiểm định ý nghĩa thống kê. Theo lý thuyết thống kê, việc kiểm định độ lớn các thông số trong mô hình thường sử dụng công thức sau [6]: ni i ni s n t  (2.18) 36 Trong đó:tni là chỉ tiêu kiểm định độ lớn của thông số thứ i, ni là giá trị của thông số thứ i, sni là độ lệch chuẩn của thông số thứ i. - Nếu tni   tα/2,n-1, (với = 0.05), thông số đó sẽ được giữ lại; - Nếu tni  < tα/2,n-1, thông số đó sẽ bị loại bỏ và quá trình tính toán để xác định giá trị của thông số sẽ được thực hiện lại theo các thông số được giữ lại. 3) Kiểm định tính độc lập giữa các thông số Quá trình phân tích trong bước kiểm định độ lớn của các thông số, đã loại đi nhiều trường hợp. Tuy nhiên, vẫn có thể có một số thông số có sự tương quan mật thiết với nhau. Để xem xét tính độc lập giữa các thông số, sẽ xem xét ma trận tương quan của các thông số. Theo [6], nếu giữa các thông số có tương quan cao (r>0,8- 0,9) thì sẽ loại bỏ một trong 2 thông số tạo nên hệ số tương quan lớn này và một mô hình có số thông số ít hơn sẽ được chọn. 2.2.5. Phần mềm thống kê SAS đối với mô hình ARIMA và ARIMAX Việc tính toán đối với mô hình ARIMA và ARIMAX là bài toán khá phức tạp, các bước kiểm nghiệm lặp lại nhiều lần. Để hỗ trợ trong quá trình tính toán đảm bảo độ chính xác, nhiều phần mềm thống kê như: STATISTICA, NCSS, SYSTAT, SASđều có chức năng tính toán các thông số trong mô hình ARIMA và ARIMAX. Trong luận văn này đã sử dụng phần mềm thống kê SAS trong việc tính toán các đặc trưng thống kê trong mô hình. SAS (Statistical Analysis System) là hệ thống phần mềm thống kê do viện nghiên cứu phần mềm thống kê của Mỹ xây dựng và đã được ở rất nhiều quốc gia trên thế giới sử dụng. Trong phần mềm SAS có đầy đủ các chức năng thống kê như lưu trữ, quản lý, kiểm tra, phân tích dử liệu. Đối với mô hình ARIMA và ARIMAX trong phần mềm SAS bao gồm một hệ thống các câu lệnh được thực hiện trên mã nguồn mở nhằm tính toán các tham số trong mô hình. Chi tiết về cách sử dụng phần mềm này được trình bày trong [37,5]. 2.3. Các nguồn số liệu đƣợc sử dụng Trên cơ sở các mục tiêu và nội dung nghiên cứu luận văn, đã sử dụng các loại số liệu sau: 1) Số liệu quan trắc lượng mưa tại 9 trạm đại diện cho khu vực 37 đồng bằng Bắc Bộ; 2) số liệu về các chỉ số khí hậu; 3) Số liệu vết đen Mặt Trời. Trong đó, chuỗi số liệu từ năm 1951 đến 2008 để phát triển mô hình, từ năm 2009 đến 2013 để kiểm chứng mô hình. Chi tiết về các nguồn số liệu và cách xử lý đối với từng loại số liệu sẽ được trình bày dưới đây: 2.3.1. Số liệu quan trắc mưa từ các trạm khí tượng Bộ số liệu tổng lượng mưa tháng được thu thập tại Trung tâm Tư liệu Khí tượng Thủy văn, hầu hết chuỗi dữ liệu được thu thập từ năm 1961 -2013, (riêng trạm Hà Nội, từ năm 1951 -2013). Nhìn chung, nguồn số liệu này có độ tin cậy cao do đã được kiểm tra, kiểm soát. Bộ số liệu này không chỉ có vai trò làm nhân tố dự báo mà còn được sử dụng làm yếu tố dự báo trong mô hình ARIMA và ARIMAX. 2.3.2. Số liệu về các chỉ số khí hậu Bộ số liệu các chỉ số khí hậu tham gia làm nhân tố dự báo bao gồm: chỉ số dao đông nam (Southern Oscillation Index, SOI), dị thường nhiệt độ mặt nước biển (ASST) trên các vùng NINO1.2, NINO3, NINO4, NINO3.4, vị trí của nhóm nhân tố này được trình bày trên hình 2.5. Các nhân tố này phản ánh khá đầy đủ hoạt động của hiện tượng ENSO trên khu vực xích đạo Thái Bình Dương. Nước ta nằm trong khu vực chịu ảnh hưởng của hiện tượng này, do vậy diễn biến của 5 đặc trương khí tượng hải dương này sẽ có quan hệ với diễn biến của lượng mưa trên các vùng lãnh thổ nước ta. Bộ số liệu này có từ năm 1951 đến nay được cập nhật thường xuyên qua mạng Internet thông qua website: edoc. Hình 2.5 [2] Vị trí nhóm nhân tố ENSO Tahiti Darwin SOI 38 Trong các chỉ số khí hậu thì nhóm nhân tố ENSO có quan hệ rất chặt chẽ với nhau và gần như chúng có chung chu kỳ (hình2.6), do vậy, tại một bước trễ thời gian, các nhân tố dự báo có thể chỉ có từ một đến 2 nhân tố thỏa mãn ý nghĩa thống kê tham gia vào mô hình dự báo. Hình 2.6. [9] Mật độ phổ của chỉ số SOI và nhiệt độ bề mặt nước biển ở các vùng Nino. 2.3.3. Số liệu về số vết đen mặt trời (Sunspot Number) Theo nhiều nghiên cứu [40] vết đen mặt trời có ảnh hưởng đến nhiều hiện tượng vật lý khí quyển, trong đó ảnh hưởng rõ rệt nhất đến lượng mưa, hạn hán và lũ lụt theo các chu kỳ và các thời gian khác nhau, số vết đen mặt trời trong một năm nào đó thường ảnh hưởng đến lượng mưa ở nhiều năm sau đó. Với ý nghĩa quan trọng như vậy, luận văn đã sử dụng số liệu này như là một trong những nhân tố dự báo lượng mưa. Nguồn số liệu vết đen mặt trời được công bố bởi Trung tâm phân tích ảnh hưởng mặt trời (Solar Influences Data analysis Center), thuộc đài quan sát hoàng gia của Bỉ,. Đây là bộ số liệu được quan trắc hàng ngày từ năm 1818 đến nay và được cập nhật thường xuyên trên trang webside: data/dailyssn.php. 2.3.4. Xử lý số liệu Trên cơ sở các nguồn số liệu được thu thập trong thời kỳ từ 1951 đến 2012, tiến hành xử lý, tính toán theo các bước sau: 39 Bước 1: Tính toán bộ số liệu về tổng lượng mưa vụ đông xuân (từ tháng 10 năm trước đến tháng 1 năm sau) đối với 9 trạm khí tượng đại diện cho khu vực đồng bằng Bắc Bộ. Đối với bộ số liệu nhân tố dự báo, nhóm các nhân tố ENSO là các đặc trưng theo tháng ở dạng chuẩn sai nên được dữ nguyên, số liệu về số vết đen Mặt trời được tính trung bình theo tháng. Bước 2: Chuẩn hóa tập nhân tố dự báo Do các nhân tố dự báo không có cùng thứ nguyên, ngoại tr