Áp dụng tính đối xứng của độ thô nháp hàm mô hình để tính toán kiểm định điểm gãy trong phân tích hồi quy - Phan Thu Hà

Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số 37, 06 - 2015 125 ÁP DỤNG TÍNH ĐỐI XỨNG CỦA ĐỘ THÔ NHÁP HÀM MÔ HÌNH ĐỂ TÍNH TOÁN KIỂM ĐỊNH ĐIỂM GÃY TRONG PHÂN TÍCH HỒI QUY Phan Thu Hà* Tóm tắt: Bài báo điểm lại một loạt phương pháp phân tích điểm gãy hiệu quả hay được sử dụng, trong đó có phương pháp tiêu chuẩn thông tin Schwartz SIC của Chen, phương pháp kiểm định sự tồn tại hệ số xu thế của Jaruskova, phương pháp tỷ số hợp lý số ELR của Liu và Qian, phương pháp tỷ số hợp lý số hiệu chỉnh MELR của Zhao, Chen và Ning cũng như một số ưu - nhược điểm của chúng. Tính đối xứng của độ thô nháp của hàm mô hình gãy khúc liên tục được thảo luận. Các phương pháp trên được áp dụng trong phân tích dữ liệu bảo hiểm. Từ khóa: Mô hình hồi quy, Điểm gãy, Tỷ số hợp lý, Độ thô nháp. 1. GIỚI THIỆU Vấn đề điểm gãy nhận được sự quan tâm ngày càng lớn trong phân tích hồi quy. Điểm gãy ngụ ý rằng, cấu trúc của dữ liệu đã có sự biến đổi kể từ thời điểm nào đó trong quá trình quan sát dài hạn. Nếu xảy ra điều này, mô hình gốc không còn phù hợp để dự báo. Khi ấy, người ta cần hiệu chỉnh tính thuần nhất của dữ liệu, dẫn đến mô hình dự báo chính xác hơn. Chẳng hạn, cần quan tâm đến xác định sớm thời điểm bắt đầu một cuộc suy thoái, hay bắt đầu một thời kỳ tăng trưởng. Xét mô hình hồi quy tuyến tính hai pha T T i i (0, k] i (k, n] iy I (i) I (i) , i 1,...,n,    x x  (1) trong đó pi x  là biến giải thích p-chiều, p,    là tham số véc tơ, 1 k n  là thời điểm chuyển chưa biết, tại đó mô hình chuyển từ trạng thái này sang trạng thái khác, i{ } độc lập, cùng phân bố, 2 i iE( ) 0, Var( )     chưa biết, AI (.) - hàm chỉ tiêu của A. Vấn đề quan tâm là đã xảy ra chuyển hay chưa, và nếu đã xảy ra thì xảy ra khi nào. Về vấn đề thứ nhất, theo ngôn ngữ của thống kê, chúng ta phải xét bài toán kiểm định: 0H :   đối với đối thuyết 1H :   . Để thuận lợi, chúng ta viết mô hình (1) dưới dạng ma trận. Muốn vậy, đặt 1 k 1 1 1 1k 1k 2 k 2 k k n n n y y y , ... , ... , ... , ... , y y y                                                      y y y y y     T T k 1 1k kp 1k 2 k k ( n k ) p 2 kT T k n ... , ... , .                               x x X O X X X = O X x x 1 trong đó stO là ma trận không cỡ s t . Lưu ý rằng khi k n thì 2k , y 2k  X và không có điểm gãy. Ta viết lại (1) như sau: k y X   (2) Với mỗi k đã cho, ước lượng bình phương cực tiểu của tham số    dựa vào đầy đủ n quan sát, k quan sát đầu và n k quan sát cuối là T 1 T T 1 T T 1 T k k k k 1 k 1k 1 k 1 k k 2 k 2 k 2 k 2 k ˆˆ ˆ) ; ) ; )  X X X y X X X y X X X y         . Các phần dư và ước lượng cho phương sai tương ứng là: T T T i i i i i i k i i i k ˆˆˆe = y x γ, e (k) = y x α , 1 i k, e (k) = y x β , k +1 i n;       n k n 2 2 2 2 2 2 i 1 i 2 i i =1 i =1 i k 1 1 1 1 ˆ ˆ ˆσ = ( e ) , σ = (e (k)) , σ = (e (k)) . n p k p n k p         Công nghệ thông tin & Khoa học máy tính Phan Thu Hà, “Áp dụng tính đối xứngtrong phân tích hồi quy.” 126 Dưới giả thuyết 0H , kˆ là tương đồng với kˆ . Tuy nhiên, dưới đối thuyết 1H , hai ước lượng này khác biệt lớn. Để tăng sức mạnh cho kiểm định, người ta đã đề nghị sử dụng phần dư dựa vào hoán đổi các tham số: T i i i k T i i i k ˆe (k) y x β , 1 i k, ˆe (k) y x α , k 1 i n.            (3) Mục 2 giới thiệu những phương pháp hay sử dụng nhất với bài toán nêu trên, đó là phương pháp tỷ số hợp lý của Quandt (1960), phương pháp SIC của Chen (1998), phương pháp phát hiện hệ số xu thế của Jarusková (1998), phương pháp tỷ số hợp lý số ELR của Liu và Quan (2010), phương pháp tỷ số hợp lý số hiệu chỉnh MELR của Zhao và ... (2013). Mục 3 đưa ra tính đối xứng của độ nháp của hàm mô hình gãy khúc liên tục có nhiều ứng dụng để tính sức mạnh của kiểm định. Mục 4 dành cho áp dụng những điều trình bày trên để phân tích dữ liệu bảo hiểm tại tỉnh V. Cuối cùng là phần kết luận. 2. MỘT SỐ PHƯƠNG PHÁP KIỂM ĐỊNH ĐIỂM GÃY a. Phương pháp tỷ số hợp lý cổ điển của Quandt ([8]) Khi các sai số  ie , i 1,...,n là độc lập, cùng phân bố chuẩn với kỳ vọng 0, độ lệch chuẩn 1 nếu i k  và 2 nếu i k  , thống kê kiểm định dựa vào tỷ số hợp lý là 3 k n 3 max { (k)}       , với k n k 1 2 n ˆ ˆ(k) (k) (k) 2 log ˆ           . (4) Giá trị lớn của  nghĩa là đã xảy ra chuyển. Tuy nhiên, Feder (1975) đã nhận xét rằng  không có phân bố 2 , Lund và Leeves (2002) chỉ ra rằng các thành phần của  (k) không độc lập. Ngày nay phương pháp này ít được sử dụng. b. Phương pháp tiêu chuẩn thông tin Schwartz SIC của Chen ([4]) Theo phương pháp này, cần phải tính các thống kê SIC(n) dưới giả thuyết 0H , và đối với k : 3 k n 2   các thống kê SIC(k) dưới đối thuyết 1H : 20 ˆ ˆSIC(n) 2 ln L ( , ) 3ln n,    2 1 0 1 0 1 ˆ ˆˆ ˆ ˆSIC(k) 2 ln L ( (k), (k), (k), (k), ) 5ln n,        (5) trong đó 20 ˆ ˆL ( , ) là hàm hợp lý cực đại ước lượng dưới giả thuyết không có chuyển ; 2 1 0 1 0 1 ˆ ˆˆ ˆ ˆL ( (k), (k), (k), (k), )     là hàm hợp lý cực đại ước lượng dưới đối thuyết có chuyển. Quy tắc quyết định là: Chọn mô hình không có chuyển nếu SIC(n) SIC(k) với mọi k; chọn mô hình với chuyển tại kˆ nếu: ˆSIC(k ) min{S(k) : 3 k n 3} SIC(n).      (6) c. Phương pháp phát hiện hệ số xu thế của Jarusková ([5]) Xét trường hợp đặc biệt của (1) khi p 2 , hàm mô hình liên tục, trong giai đoạn đầu không có hệ số xu thế, ở giai đoạn sau có thể có hệ số xu thế. Như vậy, (1) trở thành i (0, k] (k, n] iy I (i) ( b(i k) I (i) , i 1, ...,n.         (7) Với 0 k n 1   , đặt n ii 1 X (1 / n) x    ,   n 2 ii k 1 k 1/22 2 (X X) (i k) Bˆ (n k)(n k 1)(2n 2k 1) / 6 (n k) (n k 1) / 4n                , Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số 37, 06 - 2015 127 n 1 3 u (x) 2ln ln n ln x , x , 42ln ln n           s ln ( ln(1 )), t ln ( ln 1 ).          (8) Dưới những điều kiện chính quy, với n đủ lớn, giả thuyết b 0 bị bác bỏ (khi đó kết luận có điểm gãy) ở mức  nếu: k n 0 k n 1 ˆ| B | max u (t ) ˆ      (với đối thuyết hai phía | b | 0 ). (9) d. Phương pháp tỷ số hợp lý số ELR của Liu và Qian ([7]) Đối với trường hợp p 2 , thống kê n p k n p M max { 2ln (k)}       , (10) trong đó: n n i i i i i i 1 i 1 (k) sup n | e (k) 0, 0, 1                     (11) được gọi là thống kê tỷ số hợp lý số (empirical likelihood ratio), ký hiệu là ELR. Thông qua mô phỏng, Liu và Qian đã chỉ ra rằng, dưới giả thiết 0H :  n nZ = M có phân bố xấp xỉ phân bố giá trị cực trị Gumbel (là trường hợp đặc biệt của phân bố giá trị tiệm cận tổng quát) với hàm phân bố   GF (x; , ) exp exp (x ) /       , trong đó  là tham số định vị, 0  là tham số tỷ lệ; giá trị trung bình và phương sai lần lượt là C   và 2 2 / 6  , C 0.57721... là hằng số Euler. Các ông cũng chỉ ra một số giá trị ngưỡng mô phỏng của nZ với 0.10, 0.05,  0.01.  Các tham số  và 2 ứng với thống kê nZ là vững với các dạng khác nhau của phân bố của sai số. Hai tham số này tăng lên khi n tăng lên.  Dưới giả thiết 1H , phương pháp ELR có sức mạnh cao khi kích thước mẫu lớn. Để ước lượng vị trí xảy ra chuyển, phương pháp ELR tỏ ra tốt hơn phương pháp SIC của Chen, nhất là khi sự thay đổi của tham số giữa 2 pha là nhỏ hoặc khá nhỏ, và  là lớn. e) Phương pháp tỷ số hợp lý số hiệu chỉnh MELR (modified empitical likelihood ratio) của Zhao, Chen và Ning ([6]) Thống kê nM theo (10) vẫn được sử dụng. Tuy nhiên Zhao và ... đã chỉ ra rằng, dưới những điều kiện chính quy, nếu giả thuyết 0H :    là đúng thì    tn p n lim P A(ln n) M t D (ln n) exp 2e ,      (12) trong đó: A(x) 2ln x , pD (x) 2ln x (p / 2) ln ln x ln (p / 2)    . (13) Ưu điểm nổi bật của phương pháp này là dựa vào (12) cho phép tìm p-giá trị. Nhược điểm căn bản của 2 phương pháp ELR và MELR là chỉ áp dụng cho mô hình hồi quy gãy khúc liên tục, ở đó 0 1 0 0 ˆ ˆˆ ˆ k k       , nói cách khác, 0 1 0 1, , h, x f (x) x h(x x ) , x [a, b]            , (14) trong đó, các hằng số 0 1, , h, x   cho trước, h 0, x (a,b)  . Công nghệ thông tin & Khoa học máy tính Phan Thu Hà, “Áp dụng tính đối xứngtrong phân tích hồi quy.” 128 3. SỰ ĐỐI XỨNG CỦA ĐỘ THÔ NHÁP QUA TRUNG ĐIỂM Nghiên cứu độ thô nháp có vai trò quan trọng trong phân tích điểm gãy. Đối với hệ hàm cơ sở  1,x , độ thô nháp của hàm f (x), x [a,b] là giá trị:     0 1 22 0 1a, b, S (f ) min f (x) ( x) dx        . (15) Độ thô nháp của hàm gãy khúc liên tục không phụ thuộc vào hệ số 0 1,  và có thể chuyển về trường hợp    a,b 0,1 (xem, chẳng hạn [1]). Như vậy chỉ cần xét hàm   h, x 0,0,h,x f (x) f (x) h(x x ) , x 0,1       , x (0,1)  . (16) Bổ đề 3.1. Đối với hàm mô hình (16), độ nháp của nó thỏa mãn hệ thức 2 2 h,x h,1-x S (f ) S (f )  . (17) Chứng minh. Dễ thấy cực tiểu ở (15) đạt được tại 0 1,  là nghiệm của hệ 0 1 0 1 1,1 1,x 1,f x,1 x,x x,f                    (18) với 1 1 0 0 1 1,1 dx 1, 1,x x,1 x dx , 2         1 2 0 1 x,x x dx , 3    1 1 2 0 x h 1,f f (x)dx h(x x )dx (1 x ) , 2          1 2 x h x,f hx(x x )dx (1 x ) (2 x ). 6          Hệ (18) trở thành 2 0 1 2 0 1 2 h(1 x ) 3 2 h(1 x ) (2 x ).                 Giải ra ta có: 2 20 1 ˆ ˆh(1 x ) x , h(1 x ) (1 2x )           . (19) Vậy: x 12 2 2 0 1 0 1h, x 0 x ˆ ˆ ˆ ˆS (f ) (0 x) dx (h(x x ) x) dx                x 12 2 0 1 1 0 1 20 x ˆ ˆ ˆ ˆ( x) dx [(h )x (hx )] dx I I .                (20) Bây giờ xét điểm chuyển tại y 1 x   . Áp dụng (19) và (20) ta nhận được hệ số tối ưu 2 20 1 ˆ ˆh(1 y ) y , h(1 y ) (1 2y )           và độ thô nháp y 12 2 2 0 1 1 0 1 2h, y 0 y ˆ ˆ ˆ ˆS (f ) ( x) dx [(h )x (hy )] dx J J .                Ta sẽ chứng minh 1 2 2 1 0 1y ˆ ˆJ [(h )x (hy )] dx I .        Thực vậy, dùng phép đổi biến t 1 x  thì x 1 t, dx dt    , từ đó: Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số 37, 06 - 2015 129 0 2 2 1 01 y x 2 1 1 00 ˆ ˆJ [(h )(1 t) (hy )] ( dt) ˆ ˆ ˆ[( h)t h hy ] dt.                         Bởi vì 21 ˆ h h[(1 y ) (1 2y ) 1]       2 3 2 2 1 h[x (1 2(1 x )) 1] h( 2x 3x 1) ˆh(x 1) (1 2x ) ,                     2 2 1 0 2 2 2 2 0 ˆ ˆh hy h h(1 y ) (1 2y ) hy h(1 y ) y hx h(1 y ) ( 1 2y y ) hx hx ( 2 x ) ˆhx (1 2x x ) h(1 x ) x                                                nên x 2 2 1 0 10 ˆ ˆJ ( t ) dt I .       Tương tự 2 1I I , vậy 2 2 h, x h, y S (f ) S (f ).  Mệnh đề được chứng minh. Định lý 3.2. Độ thô nháp của hàm mô hình (14) là đối xứng qua trung điểm của đoạn quan sát [a, b] và đạt được giá trị cực đại 2h /192 tại trung điểm của đoạn này. Chứng minh. Áp dụng Bổ đề 3.1 với lưu ý rằng độ nháp cho trường hợp đoạn quan sát [a,b] bất kỳ có thể chuyển về trường hợp đoạn quan sát trên đoạn [0,1], chúng ta nhận được phần đầu của Định lý. Dùng đổi biến đưa [a, b] về [0, 1], phần sau được suy ra từ Định lý 2 trong [3]. Xét (2) với p 2 , ix i / n , hàm mô hình liên tục, k  đã biết, 2i i.i.d.N(0, }  , 2 chưa biết. Đặt T0 1ˆ ˆ ˆ( , ) ,   T 0 1k ˆ ˆ ˆ, ) ,    T 0 1k ˆ ˆ ˆ, ) .    Để kiểm định giả thuyết 0H không có điểm gãy đối với đối thuyết 1H đã xảy ra điểm gãy, có thể dùng quy tắc bác bỏ giả thuyết khi 1 2 2F (F F ) / F f (p, n 2p)    (21) trong đó:     2k n2 1 0 1 i 0 1 i 0 1 i 0 1 ii 1 i k 1 ˆ ˆˆ ˆ ˆ ˆ ˆ ˆF t ( t ) t ( t )                    ,     2k n2 2 i 0 1 i i 0 1 ii 1 i k 1 ˆ ˆˆ ˆF y ( t ) y ( t )              , f (p, n 2p)  là phân vị mức của phân bố F với p và n 2p bậc tự do. Từ Định lý 3 trong [2], sức mạnh của (21) là P{F(p,n 2p,nR) f (p,n 2p)}   , với F(p,n 2p,nR) là biến ngẫu nhiên có phân bố F phi trung tâm với p, n 2p bậc tự do và tham số phi trung tâm nR, R là độ nháp của hàm * 1 1 * 2 2 A B x, 0 x k / n, f (x) A B x, k / n x 1.          Từ Định lý 2 chúng ta nhận được: Hệ quả 3.3. Sức mạnh của kiểm định (21) là hàm đối xứng của điểm chuyển x k / n  qua trung điểm 1 / 2 của đoạn [0,1]. Nhận xét. Kết quả ở Hệ quả 3.3 đã được Farley, J.U. và cộng sự (1985) và [5] nhắc đến, song họ chỉ nêu lên hiện tượng thông qua mô phỏng. Công nghệ thông tin & Khoa học máy tính Phan Thu Hà, “Áp dụng tính đối xứngtrong phân tích hồi quy.” 130 4. NGHIÊN CỨU THỰC NGHIỆM Chúng ta sẽ dùng các phân tích trên trên để xử lý dữ liệu chi trả lương hưu và tử tuất (tỷ VNĐ) hàng tháng, từ 1/2008 đến 12/2014 tại tỉnh V ở hình 1. Như vậy, có 83 tháng, ta đánh số từ tháng thứ nhất (1/2008) đến tháng thứ 83 (12/2014). Trước hết chúng ta cần khảo sát dữ liệu ngoại lai (quá thô). Mô hình hồi quy tuyến tính đơn được dùng để lọc mô hình. Dữ liệu số 69 (9/2013) nhận giá trị 27.336 có phần dư chuẩn hóa 69d 4.02794 3  . Ta loại đi dữ liệu này. Dữ liệu khuyết còn lại thể hiện mô hình hai pha và không liên tục khá rõ. Chúng ta sẽ sử dụng phương pháp SIC của Chen để phát hiện điểm chuyển. Viết lại (6) dưới dạng       k n n2 2 22 i i i1 k 1 13 k n 3 ˆmin e (k) e (k) 2 ln n e .              (22) Ở đây 2ˆn 83, 1.6316   nên 2ˆ2 ln n 14.4196  . Cực tiểu ở vế trái (22) đạt được tại k 52  . Vế trái nhận giá trị 43.3095, nhỏ hơn rất nhiều so với vế phải 132.1629. Chúng ta kết luận có chuyển tại k 52   4/2012. Hàm hồi quy ước lượng cho giai đoạn trước và sau chuyển lần lượt là (xem hình 1): y 7.4257 0.1477t, t 52, y 8.4324 0.1984t, 52 t 84.         (23) Hình 1. Chi trả lương hưu-tử tuất, dữ liệu 69 là ngoại lai, đường đứt đứng thể hiện điểm gãy, các đường đứt xiên thể hiện xu thế trước và sau chuyển. Độ nháp của hàm này theo (15) là 0.12129, một giá trị khá lớn. Kiểm định T cho kết luận các hệ số ước lượng được là khác 0 có ý nghĩa. Cần kiểm định phương sai hai giai đoạn bằng nhau: 1 2   . Ước lượng cho sai số chuẩn tương ứng là 1ˆ 0.58129,  2ˆ 0.71907  . Từ đó 2 2 2 1ˆ ˆ/ 1.53023   0.05f (29, 50) 1.69  . Vậy chúng ta không bác bỏ giả thuyết 1 2   , nói cách khác, ta chấp nhận phương sai hai giai đoạn bằng nhau. Hình 2. Đồ thị các phần dư, đường đứt đứng chỉ điểm chuyển tại k 52  . 5. KẾT LUẬN Cần kiểm tra tính độc lập của các phần dư (Hình 2). Hàm tự tương quan của chúng là Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số 37, 06 - 2015 131 ar(1) 0 , ar(2) 0 ,.53395 ar(3).39 0182 .19593   , đều là những giá trị lớn có ý nghĩa. Như vậy, có tương quan chuỗi cao giữa các phần dư. Tóm lại, dữ liệu đã có một quan sát ngoại lai vào tháng 9/2013. Đã xảy ra chuyển tại k 52  (ứng với thời điểm 4/2012) theo mô hình (23). Ở giai đoạn sau, hệ số chặn và hệ số góc đều tăng lên đáng kể, thể hiện chi trả lương hưu-tử tuất ở mức độ hoàn toàn khác với giai đoạn trước. Khi cần dự báo cho giai đoạn sau, có thể phải dùng mô hình AR(3). Bài báo đã đưa ra và phân tích các ưu - nhược điểm của các phương pháp điểm chuyển hiệu quả đang được áp dụng. Tính đối xứng của độ nháp của hàm mô hình gãy liên tục qua trung điểm của đoạn quan sát được chứng minh, từ đó nhận được tính đối xứng của sức mạnh của kiểm định điểm gãy. Phân tích trên được áp dụng cho dữ liệu bảo hiểm. Kinh nghiệm chỉ ra rằng, điều quan trọng hàng đầu là lựa chọn dạng mô hình phù hợp. TÀI LIỆU THAM KHẢO [1] Phan Thu Hà, “Mức độ thô nháp của hàm số và ứng dụng”, Tạp chí Nghiên cứu KH & CN Quân sự, 30, (2014), pag.34-39. [2] Ban, T.V., Quyen, N.T., The power of change-point test for two-phase regression, Applied Mathematics, 5, (2014), pag.2994-3000. [3] Ban, T.V., Quyen, N.T., Ha, P.T., “The roughness of model function to the basis functions”, J. Math. and System Science, 3, No. 8 (2013), pag.385-390. [4] Chen, J., “Testing for a change point in linear regression models”, Comm. in Statistics.-Theory and Methods, 27, No.10 (1998), pag.2481-2493. [5] Jarusková, J., “Testing appearance of linear trend”, J. Statistical Planning and Inference, 70 (1998), pag.263 – 276. [6] Zhao, H., Chen, H., and Ning, W., “Changepoint Analysis by Modified Empirical Likelihood Method in Two-phase Linear Regression Models”, Open Journal of Applied Sciences, 3 (2013), pag.1-6. [7] Z. Liu and L. Qian, “Changepoint Estimation in a Segmented Linear Regression via Empirical Likelihood”, Communications in Statistics-Simulation and Computation, 89 (2010), pag.85-100. [8] Quandt, R. E., “Tests of the hypothesis that a linear regression system obeys two separate regimes”, J. Amer. Stat. Ass. 55 (1960), pag.324–330. ABSTRACT APPLYING THE SYMMETRY OF ROUGHNESS OF THE MODEL FUNCTION FOR CALCULATING CHANGE-POINT TESTS IN REGRESSION The article points out the series of efficient change-point analytical methods including the Chen’s Schwartz information criteria (SIC) method, the Jaruskova’s methods above testing appearance of linear trend, the Liu & Qian’s empirical likelihood (ELR) method, the Zhao, Chen & Ning’s modified empirical likelihood method (ELR) as well as a number of advantages – disadvantages of them. Symmetry of the roughness of the continuous segment model function is discussed. These methods are applied in the analysis of insurance data. Keywords: Regression model, Changepoint, Likelihood ratio, Roughness. Nhận bài ngày 04 tháng 3 năm 2015 Hoàn thiện ngày 02 tháng 4 năm 2015 Chấp nhận đăng ngày 12 tháng 06 năm 2015 Địa chỉ: Khoa Công nghệ Thông tin, Học viện KTQS, ĐT: 0985 193 986 *Email: phanthuha_bmt@yahoo.com