Phân tích chuỗi thời gian với sự hỗ trợ của package analyzets

ndex da ta 0 10 20 30 40 10 0 15 0 20 0 ts forecast Chen-Hsu 10 fuzzy set index da ta 0 10 20 30 40 10 0 15 0 20 0 ts forecast Chen-Hsu 12 fuzzy set index da ta 0 10 20 30 40 10 0 15 0 20 0 ts forecast Chen-Hsu 15 fuzzy set index da ta 0 10 20 30 40 10 0 15 0 20 0 ts forecast Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 44 > goc<-data.frame(ban) > chen.hsu.density<-data.frame(chen.hsu8.density, +chen.hsu10.density,chen.hsu12.density,chen.hsu15.density) > av.res(Y=goc,F=chen.hsu.density) chen.hsu8.density chen.hsu10.density chen.hsu12.density ME 1.431 0.065 0.586 MAE 5.09 3.635 3.75 MPE 0.9 -0.109 0.478 MAPE 3.989 2.761 2.941 MSE 42.539 21.814 20.911 RMSE 6.522 4.671 4.573 U 0.093 0.066 0.065 chen.hsu15.density min.model ME 0.227 chen.hsu10.density MAE 3.069 chen.hsu15.density MPE 0.062 chen.hsu10.density MAPE 2.389 chen.hsu15.density MSE 13.764 chen.hsu15.density RMSE 3.71 chen.hsu15.density U 0.053 chen.hsu15.density Theo tiêu chuẩn MAPE mô hình Chen-Hsu với 15 tập mờ (chen.hsu15.density) là tốt nhất. Mô hình Chen-Hsu với bin != NULL Nếu bin = NULL thì số tập mờ sẽ được R chia lại lần hai theo giá trị của tham số‘divide’. Khi đó, số lượng phần tử có trong các tập mờ iA , ni ,1 có thể không như ý muốn. Để R chia lại số tập mờ theo ý muốn, ta dùng tham số ‘bin’ của hàm. Cho ‘bin’ một vector số mà các phần tử của vector là các điểm chia của các tập mờ (không có giá trị của điểm đầu và điểm cuối) khi đó R sẽ chia số tập mờ theo tham số ‘bin’. Trong lần phân tập mờ đầu tiên, chúng tôi cho số tập mờ là 8 và được kết quả của 8 tập mờ. > fuzzy.ts1(ban,n=8,type="Chen-Hsu",trace=1)$table1 set dow up mid num 1 A1 61.94 84.30625 73.12313 5 2 A2 84.30625 106.6725 95.48937 5 3 A3 106.6725 129.0388 117.8556 7 4 A4 129.0388 151.405 140.2219 8 5 A5 151.405 173.7713 162.5881 9 6 A6 173.7713 196.1375 184.9544 5 7 A7 196.1375 218.5038 207.3206 5 8 A8 218.5038 240.87 229.6869 3 Từ bảng kết quả ta thấy, các phần tử trong tập mờ thứ 4 và 5 nhiều hơn số Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 45 phần tử của các tập mờ còn lại rất nhiều. Do đó, ta sẽ chia như sau > chen.hsu.bin<-fuzzy.ts1(ban,type="Chen-Hsu",bin=c(84.30625,106.6725,129.0388, + 140.2219,151.405,162.5882,173.7713,196.1375,218.5038),plot=TRUE) Chen-Hsu 5 fuzzy set index d a ta 0 10 20 30 40 1 0 0 1 5 0 2 0 0 ts forecast Đoạn Số phần tử trong đoạn A1 = [61.94; 84.30625] 5 A2 = [84.30625; 106.6725] 5 A3 = [106.6725; 129.0388] 7 A4.1 = [129.0388;140.2219] 3 A4.2 = [140.2219;151.405] 5 A5.1 = [151.405; 162.5882] 5 A5.2 = [162.5882; 173.7713] 4 A6 = [173.7713; 196.1375] 5 A7 = [196.1375; 218.5038] 5 A8 = [218.5038; 240.87] 3 Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 46 e) Tìm mô hình mờ tốt nhất và xây dựng mô hình ARIMA > goc<-data.frame(ban) > mo<-data.frame(chen15,singh15,heuristic15,chen.hsu15.dt, + chen.hsu15.density,chen.hsu.bin) > av.res(Y=goc,F=mo) chen15 singh15 heuristic15 chen.hsu15.dt ME 1.54 -0.126 0.235 0.235 MAE 27.27 3.033 23.24 3.04 MPE -6.067 -0.206 -5.049 0.068 MAPE 22.048 2.365 18.091 2.37 MSE 1193.914 13.928 943.975 13.661 RMSE 34.553 3.732 30.724 3.696 U 0.491 0.053 0.437 0.053 chen.hsu15.density chen.hsu.bin min.model ME 0.227 1.831 singh15 MAE 3.069 4.868 singh15 MPE 0.062 1.191 chen15 MAPE 2.389 3.829 singh15 MSE 13.764 37.641 chen.hsu15.dt RMSE 3.71 6.135 chen.hsu15.dt U 0.053 0.087 chen.hsu15.dt Theo tiêu chuẩn MAPE mô hình Singh với 15 tập mờ (singh15) là tốt nhất trong tất cả các mô hình mờ. Chúng ta chọn chuỗi này để đi xây dựng mô hình ARIMA. Các bước xây dựng mô hình ARIMA tương tự mục 2.3.3 sau khi đã tìm được chuỗi làm trơn tốt nhất. Sau khi thực hiện qua các bước 1 và 2 ta tìm được mô hình tối ưu là SARIMA(1,1,0)x(0,0,1)[5]. Tuy nhiên khi thực hiện bước 3 để kiểm định, ta lại thấy chuỗi phần dư không là ngẫu nhiên trắng. Do đó mô hình này không thích hợp để dự báo cho tương lai. 2.4. Mô hình chuỗi thời gian mờ Abbasov-Mamedova > ban<-ts(dat$sales) > w2<-fuzzy.ts2(ban,n=5,r=7,w=2,C=0.00001,trace=TRUE,forecast=15,plot=TRUE) Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 47 > w2 $type [1] "Abbasov-Manedova" $table1 U low up Bw 1 u1 -140.76 -81.026 -110.893 2 u2 -81.026 -21.292 -51.159 3 u3 -21.292 38.442 8.575 4 u4 38.442 98.176 68.309 5 u5 98.176 157.91 128.043 $table2 point ts diff.ts 1 1 199.95 NA 2 2 195.74 -4.21 Abbasov-Mamedova, C = 1e-05 w = 2 n = 5 index d a ta 0 10 20 30 40 1 0 0 1 5 0 2 0 0 ts forecast Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 48 3 3 102.68 -93.06 45 45 82.96 8.83 46 46 240.87 157.91 47 47 151.52 -89.35 $table3 [1] NA [2] "A[2]={(0.9999989/u1),(0.9999998/u2),(1/u3),(0.9999995/u4),(0.9999983/u5)}" [3] "A[3]={(1/u1),(0.9999998/u2),(0.999999/u3),(0.9999974/u4),(0.9999951/u5)}" [45] "A[45]={(0.9999986/u1),(0.9999996/u2),(1/u3),(0.9999996/u4),(0.9999986/u5)}" [46] "A[46]={(0.9999928/u1),(0.9999956/u2),(0.9999978/u3),(0.9999992/u4),(0.9999999/u5)}" [47] "A[47]={(1/u1),(0.9999999/u2),(0.999999/u3),(0.9999975/u4),(0.9999953/u5)}" $table4 point interpolate diff.interpolate 1 4 111.25488 8.574882 2 5 171.45494 8.574937 3 6 110.33498 8.574979 42 45 82.7049 8.5749 43 46 91.53489 8.574894 44 47 249.44518 8.575181 $table5 point forecast diff.forecast 1 48 258.0203 8.575082 2 49 266.5952 8.574896 3 50 275.1702 8.575 4 51 283.7452 8.575 5 52 292.3202 8.575 $table6 [1] "A[48]={(0.9999986/u1),(0.9999996/u2),(1/u3),(0.9999996/u4),(0.9999986/u5)}" [2] "A[49]={(0.9999986/u1),(0.9999996/u2),(1/u3),(0.9999996/u4),(0.9999986/u5)}" [3] "A[50]={(0.9999986/u1),(0.9999996/u2),(1/u3),(0.9999996/u4),(0.9999986/u5)}" [4] "A[51]={(0.9999986/u1),(0.9999996/u2),(1/u3),(0.9999996/u4),(0.9999986/u5)}" [5] "A[52]={(0.9999986/u1),(0.9999996/u2),(1/u3),(0.9999996/u4),(0.9999986/u5)}" $accuracy ME MAE MPE MAPE MSE RMSE U Abbasov.Mamedova -7.465 59.492 -19.668 48.1 5031.496 70.933 1.005 > w3<-fuzzy.ts2(ban,n=5,r=7,w=3,C=0.00001,trace=TRUE,forecast=5) > w4<-fuzzy.ts2(ban,n=5,r=7,w=4,C=0.00001,trace=TRUE,forecast=5) Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 49 > w5<-fuzzy.ts2(ban,n=5,r=7,w=5,C=0.00001,trace=TRUE,forecast=5) > w6<-fuzzy.ts2(ban,n=5,r=7,w=6,C=0.00001,trace=TRUE,forecast=5) > sosanh<-rbind(w2$accuracy,w3$accuracy,w4$accuracy,w5$accuracy, + w6$accuracy) > dimnames(sosanh)[[1]]<-c("w2","w3","w4","w5","w6") > sosanh ME MAE MPE MAPE MSE RMSE U w2 -7.465 59.492 -19.668 48.1 5031.496 70.933 1.005 w3 -8.839 59.675 -20.862 48.482 5086.53 71.32 1.008 w4 -7.39 59.437 -19.728 48.005 5091.986 71.358 1.005 w5 -9.439 59.018 -21.209 48.177 5073.051 71.225 1.009 w6 -7.791 58.61 -20.172 47.814 5057.993 71.12 1.006 Theo tiêu chuẩn MAPE chúng ta thấy mô hình w6 (mô hình với tham số w = 6) có giá trị nhỏ nhất, ta nói mô hình Abbasov-Mamedova với tham số w=6 là đáng tin cậy nhất. Truy xuất giá trị dự báo và vẽ đồ thị dự báo của mô hình như sau: > dubao.abbasov<- w6$table5 point forecast diff.forecast 1 48 258.0201 8.57486 2 49 266.5951 8.575 3 50 275.1701 8.575 4 51 283.7451 8.575 5 52 292.3201 8.575 > plot(dubao.abbasov [,2],type="o",pch=18,col=2,xlab="point",ylab="sales", + main="Đồ thị dự báo bằng mô hình Abbasov-Mamedova") > grid.on(v=FALSE) 1 2 3 4 5 2 6 0 2 7 5 2 9 0 Ðồ thị dự báo bằng mô hình Abbasov-Mamedova point s a le s Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 50 2.5. Tổng hợp và so sánh các mô hình dự báo Sau khi phân tích qua các mô hình, ta có được 3 mô hình có thể dự báo được là mô hình với biến giả nhiệt độ, mô hình với số liệu làm trơn và mô hình Abbasov- Mamedova. Ta có kết quả dự báo từ 3 mô hình này lại thành một data frame để tiện quan sát. > kqth<-data.frame(bien.gia=dubao.bgia$pred,tron=dubao.tron$pred, + chuoi.mo=dubao.abbasov$forecast) > kqth bien.gia tron chuoi.mo 1 120.6013 150.7967 258.0201 2 146.9638 132.7841 266.5951 3 138.5314 122.5042 275.1701 4 143.7891 138.9585 283.7451 5 139.8764 156.1941 292.3201 Nhìn chung ngoài mô hình Abbasov-Mamedova có giá trị dự báo tăng liên tục, thì 2 mô hình còn lại đều có giá trị tăng giảm không đều trong 5 ngày tiếp theo. > bien.gia<-fitted(fit.bgia) > lam.tron<-c(NA,NA,fitted(fit.tron)) > abbasov<-c(rep(NA,7),w6$table4$interpolate) > noisuy<-data.frame(bien.gia,lam.tron,abbasov) > goc<-data.frame(ban) > av.res(Y=goc,F=noisuy) bien.gia lam.tron abbasov min.model ME -3.628 -1.087 -7.791 abbasov MAE 33.414 37.92 58.61 bien.gia MPE -12.893 -11.267 -20.172 abbasov MAPE 28.259 30.344 47.814 bien.gia MSE 1771.301 2061.165 5057.993 bien.gia RMSE 42.087 45.4 71.12 bien.gia U 0.605 0.638 1.006 bien.gia Theo tiêu chuẩn MAPE thì mô hình bien.gia (mô hình ARIMA với biến giả nhiệt độ) là đáng tin cậy nhất. Các lệnh bên dưới giúp vẽ một đồ thị tổng hợp cho kết quả dự báo. > plot(kqth$bien.gia,type="o",pch=16,col="blue",xlab="point",ylab="sales", + main="Kết quả dự báo từ các mô hình",ylim=c(min(kqth),max(kqth))) > lines(kqth$tron~c(48:52),pch=17,type="o",col="green") > lines(kqth$chuoi.mo~c(48:52),pch=18,type="o",col="red") > legend(51,250,c("biến giả","làm trơn","chuỗi mờ"),lty=c(1,1,1), + pch=c(16,17,18),col=c("blue","green","red")) > grid.on() Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 51 3. Bài toán tỷ xuất sinh lợi giá cổ phiếu 3.1. Nguồn số liệu Trong phân tích tài chính, nếu kí hiệu một loại giá cả thay đổi hàng ngày là ,tGia để thuận tiện cho việc phân tích các yếu tố ngẫu nhiên của các chỉ số, người ta hay xét các đại lượng 1 1 t t t t Gia Gia SST Gia     Hay 1ln lnt t tSSL Gia Gia   được xem là tỷ suất sinh lợi (lợi nhuận) hay lợi nhuận logarit (t 1) , tGia là giá tại thời điểm t. Giả thiết thường được ưa chuộng nhất là giả thiết cho rằng 1,..., tSSL SSL tuân theo luật phân phối chuẩn. Thế nhưng, theo sự phân tích các chuỗi thời gian tài chính, giả thiết đó rất nhiều khi không phù hợp với thực tế diễn biến của giá cả tài chính. Ta có thể phân tích chuỗi giá chứng khoán thông qua việc phân tích chuỗi .tSSL Dữ liệu phân tích cho bài toán này là giá cổ phiếu ACB – mã cổ phiếu của Ngân hàng Á Châu. Nguồn dữ liệu được lấy từ thời điểm 21/11/2006 đến 29/09/2015 tại địa chỉ www.Cophieu68.com/export/excel.php?id=ACB Kết quả dự báo từ các mô hình point sa le s 48 49 50 51 52 1 5 0 2 0 0 2 5 0 biến giả làm tron chuỗi mờ Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 52 Giá cổ phiếu được tính lúc đóng cửa (đơn vị ngàn đồng). 3.1. Mô hình ARMA–GARCH Bước 1: Thiết lập suất sinh lợi và xác định mô hình dự báo Ta tạo ra chuỗi dữ liệu suất sinh lợi với 2199 quan sát, trong đó suất sinh lợi của chỉ số giá chứng khoán ACB được tính theo công thức SSL(t)=ln(Gia(t)/Gia(t-1)). > data<-read.csv(file.choose(),header=TRUE) > names(data) [1] "X.Ticker." "X.DTYYYYMMDD." "X.OpenFixed." "X.HighFixed." [5] "X.LowFixed." "X.CloseFixed." "X.Volume." "X.Open." [9] "X.High." "X.Low." "X.Close." "X.VolumeDeal." [13] "X.VolumeFB." "X.VolumeFS." "date" "month" > head(data1) X.DTYYYYMMDD. X.CloseFixed. date month 1 20150929 19.4 29 9 2 20150928 19.3 28 9 3 20150925 19.7 25 9 4 20150924 19.5 24 9 5 20150923 19.3 23 9 6 20150922 19.2 22 9 Đây là file số liệu gốc, nên bộ số liệu có khá nhiều cột và thứ tự các ngày không đúng với thứ tự chuỗi thời gian. Nên ta cần đảo lại thứ tự của dữ liệu và trích lọc chuỗi giá đóng cửa (cột X.CloseFixed) để phân tích. Đoạn code phía sau làm 2 công việc này. > temp<-data1 > data2<-temp > for(i in 1:dim(data2)[1]) data2[i,]<-temp[length(data[,6])-i+1,] > head(data2) X.DTYYYYMMDD. X.CloseFixed. date month 1 20061121 26.3 21 11 2 20061122 27.2 22 11 3 20061123 27.9 23 11 4 20061124 29.6 24 11 5 20061127 29.6 27 11 6 20061128 28.9 28 11 > close<-ts(data2[,2],start=1,frequency=5) Chuỗi giá đóng của đã được lưu trữ trong đối tượng close dưới dạng chuỗi thời gian, ta có thể dùng nó để phân tích. > #Suất sinh lợi > ln.close<-log(close) > SSL<-diff(ln.close) > Descriptives(SSL,plot="TRUE",r=0,answer=2) Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 53 N: NaN: Min: 1sq QU: Median: Mean: 3rd QU: Max: VAR: SD: SE: x 2198 0 0 0 0 0 0 0 0 0 0 > library(urca) > summary(ur.df(SSL,type="none",lags=5,selectlag="BIC")) Tóm tắt kết quả kiểm định ADF. Test value Giá trị tới hạn tau 1% 5% 10% -32.0015 -2.58 -1.95 -1.62 Chuỗi SSL dao động quanh giá trị 0, kiểm định ADF cho giá trị thống kê |tau| = |-32.0015| > |-1.95| (giá trị tới hạn tau 5%). Ta có đủ bằng chứng để kết luận chuỗi SSL là chuỗi dừng. Từ biểu đồ ACF và PACF chúng ta xác định được mô hình cần tìm là SARIMA với các bậc là p = 2, d = 0, q = 1, P = 1, D = 0, Q = 1. Time x 0 100 200 300 400 -0 .1 0 0 .0 0 0 .1 0 Histogram of x x D e n s it y -0.10 0.00 0.10 0 1 0 2 0 3 0 4 0 5 0 Skewness: 0 Kurtosis: 4 0 1 2 3 4 5 6 0 .0 0 0 .0 5 0 .1 0 0 .1 5 Series x Lag A C F 0 1 2 3 4 5 6 -0 .0 5 0 .0 5 0 .1 5 Lag P a rt ia l A C F Series x Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 54 Bước 2: Tìm mô hình tối ưu > PrintAIC(SSL,order=c(2,0,1),seas=list(order=c(1,0,1),frequency=5),type="SARIMA") $mohinh Mo hinh Gia tri AIC Xep loai mo hinh 1 SARIMA(0,0,1)*(0,0,1),s=5 ...AIC = -10306.5 .......... 1 mo hinh 2 SARIMA(0,0,1)*(1,0,0),s=5 ...AIC = -10305.59 .......... 4 mo hinh 3 SARIMA(0,0,1)*(1,0,1),s=5 ...AIC = -10306.31 .......... 2 mo hinh 4 SARIMA(1,0,0)*(0,0,1),s=5 ...AIC = -10304.28 .......... 9 mo hinh 5 SARIMA(1,0,0)*(1,0,0),s=5 ...AIC = -10303.32 ......... 15 mo hinh 6 SARIMA(1,0,0)*(1,0,1),s=5 ...AIC = -10304.18 ......... 11 mo hinh 7 SARIMA(1,0,1)*(0,0,1),s=5 ...AIC = -10304.54 .......... 7 mo hinh 8 SARIMA(1,0,1)*(1,0,0),s=5 ...AIC = -10303.64 ......... 13 mo hinh 9 SARIMA(1,0,1)*(1,0,1),s=5 ...AIC = -10304.35 .......... 8 mo hinh 10 SARIMA(2,0,0)*(0,0,1),s=5 ...AIC = -10305.6 .......... 3 mo hinh 11 SARIMA(2,0,0)*(1,0,0),s=5 ...AIC = -10304.74 .......... 6 mo hinh 12 SARIMA(2,0,0)*(1,0,1),s=5 ...AIC = -10305.35 .......... 5 mo hinh 13 SARIMA(2,0,1)*(0,0,1),s=5 ...AIC = -10304.18 ......... 10 mo hinh 14 SARIMA(2,0,1)*(1,0,0),s=5 ...AIC = -10303.34 ......... 14 mo hinh 15 SARIMA(2,0,1)*(1,0,1),s=5 ...AIC = -10303.89 ......... 12 $best Mo hinh toi uu SARIMA(0,0,1)*(0,0,1),s=5 ...AIC = -10306.5 Mô hình tối ưu là SARIMA(0,0,1)*(0,0,1)[5]. Ước lượng hệ số cho mô hìnhtối ưu. > fit1<-arima(SSL,order=c(0,0,1),seas=list(order=c(0,0,1),5)) > fit1 Call: arima(x = SSL, order = c(0, 0, 1), seasonal = list(order = c(0, 0, 1), 5)) Coefficients: ma1 sma1 intercept 0.1514 0.0903 -1.00E-04 s.e. 0.0211 0.0219 6.00E-04 sigma^2 estimated as 0.0005364: log likelihood = 5157.25, aic = -10308.5 Bước 3: Kiểm tra hiệu ứng ARCH > Descriptives(res,answer=2,plot=TRUE) N: NaN: Min: 1sq QU: Median: Mean: 3rd QU: Max: VAR: SD: SE: x 2198 0 -0.12 -0.01 0 0 0.01 0.13 0 0.02 0 Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 55 > Box.test(res,lag=32,type="Ljung") Box-Ljung test data: res X-squared = 47.459, df = 32, p-value = 0.03855 Chuỗi phần dư có các giá trị dao động quanh giá trị 0, nhưng có phương sai không bằng nhau. Trên biểu đồ ACF có nhiều độ trễ lớn hơn 0, kiểm định Ljung – Box cho giá trị p–value = 1.00e-04< 0.05, cho thấy có hiện tượng tương quan chuỗi trong chuỗi phần dư. Ta thực hiện kiểm định nhân tử Lagrange để xem có ảnh hưởng ARCH/GARCH hay không. Time x 0 100 200 300 400 -0 .1 0 0 .0 0 0 .1 0 Histogram of x x D e n s it y -0.10 0.00 0.10 0 1 0 2 0 3 0 4 0 Skewness: 0.15 Kurtosis: 3.74 0 1 2 3 4 5 6 -0 .0 4 0 .0 0 0 .0 4 Series x Lag A C F 0 1 2 3 4 5 6 -0 .0 4 0 .0 0 0 .0 4 Lag P a rt ia l A C F Series x Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 56 > library(FinTS) > ArchTest(res,30) ARCH LM-test; Null hypothesis: no ARCH effects data: res Chi-squared = 418.5, df = 30, p-value < 2.2e-16 Giá trị p-value < 2.2e-16 (rất nhỏ hơn so với mức ý nghĩa 5%) chúng ta bác bỏ giả thuyết, tức là chuỗi phần dư có ảnh hưởng ARCH/GARCH. Bước 4: Xác định bậc của ARCH > tsnew<-res^2 > par(mfrow=c(2,1)) > acf(tsnew,lag=100) > pacf(tsnew,lag=100) 0 5 10 15 20 -0 .0 5 0 .1 5 0 .3 0 Series tsnew Lag A C F 0 5 10 15 20 -0 .0 5 0 .1 0 0 .2 5 Lag P a rt ia l A C F Series tsnew Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 57 Từ biểu đồ trên ta thấy, biểu đồ ACF tắt dần về 0 nên ta khẳng định chuỗi phần dư có ảnh hưởng ARCH. Chọn bậc q = 9, ta tìm mô hình tối ưu từ các mô hình dự tuyển. > PrintAIC(res,order=9,type="ARCH") $mohinh Mo hinh Gia tri AIC Xep loai mo hinh 1 ARCH(1) ...AIC = -10845.84 .......... 9 mo hinh 2 ARCH(2) ...AIC = -11186.86 .......... 8 mo hinh 3 ARCH(3) ...AIC = -11230.91 .......... 7 mo hinh 4 ARCH(4) ...AIC = -11345.73 .......... 6 mo hinh 5 ARCH(5) ...AIC = -11386.54 .......... 5 mo hinh 6 ARCH(6) ...AIC = -11402.35 .......... 4 mo hinh 7 ARCH(7) ...AIC = -11406.84 .......... 3 mo hinh 8 ARCH(8) ...AIC = -11451.1 .......... 2 mo hinh 9 ARCH(9) ...AIC = -11451.65 .......... 1 $best Mo hinh toi uu ARCH(9) ...AIC = -11451.65 Warning messages: 1: In sqrt(pred$e) : NaNs produced 2: In sqrt(pred$e) : NaNs produced 3: In sqrt(pred$e) : NaNs produced 4: In sqrt(pred$e) : NaNs produced 5: In sqrt(pred$e) : NaNs produced 6: In sqrt(pred$e) : NaNs produced 7: In sqrt(pred$e) : NaNs produced 8: In sqrt(pred$e) : NaNs produced Mô hình tối ưu là ARCH(8). Tuy nhiên bậc 8 là khá lớn nên ta xét thêm mô hình của GARCH(1,1) để xem chỉ số AIC của mô hình có tốt hơn không. > PrintAIC(res,order=c(1,1),type="GARCH") $mohinh Mo hinh Gia tri AIC Xep loai mo hinh 1 GARCH(1,0) ...AIC = -10306.93 .......... 2 mo hinh 2 GARCH(1,1) ...AIC = -11478.51 .......... 1 $best Mo hinh toi uu GARCH(1,1) ...AIC = -11478.51 Warning messages: Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 58 1: In garch(DataTimeSeries, order = c(i1, i2), trace = FALSE) : singular information 2: In sqrt(pred$e) : NaNs produced Chỉ số AIC của GARCH(1,1) nhở hơn só với mô hình ARCH(8). Ta chọn mô hình GARCH(1,1) để dự báo cho phương sai. Bước 5: Ước lượng mô hình Ta có bảng kết quả ước lượng của mô hình như sau: > fit2<-garch(res,order=c(1,1),trace=0) Warning message: In sqrt(pred$e) : NaNs produced > summary(fit2) Call: garch(x = res, order = c(1, 1), trace = 0) Model: GARCH(1,1) Residuals: Min 1Q Median 3Q Max -9.67232 -0.55306 -0.00464 0.538362 6.086255 Coefficient(s): Estimate Std. Error t value Pr(>|t|) a0 1.32E-05 6.76E-07 19.55 <2e-16 *** a1 2.69E-01 2.00E-02 13.46 <2e-16 *** b1 7.43E-01 1.28E-02 58.01 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Diagnostic Tests: Jarque Bera Test data: Residuals X-squared = 7449.6, df = 2, p-value < 2.2e-16 Box-Ljung test data: Squared.Residuals X-squared = 0.093516, df = 1, p-value = 0.7598 Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 59 Tất cả các hệ số của mô hình đều có ý nghĩa thống kê. Giá trị p-value của kiểm định Ljung-Box là 0.7598 > 0.05 ta nói chuỗi sai số của mô hình là nhiễu trắng. Mô hình là đầy đủ và có thể dùng để dự báo. Bước 6: Dự báo > forecastGARCH(fitGARCH=fit2,fitARMA=fit1,trace=TRUE,r=7) $ARCH a0 a1 b1 1.32e-05 2.69e-01 7.43e-01 $ARMA ma1 sma1 intercept 0.151356 0.090348 -0.00014 $forecast Point res res^2 SSL.forecast VAR.forecast 1 (440,4) 0.00892 7.96E-05 0.0002482 2 (440,5) 0.0017134 0.0002189 > sqrt( 0.0002189) [1] 0.01479527 Theo kết quả dự báo vào ngày 30/09/2015 (ngày tiếp theo của chuỗi) suất sinh lợi kỳ vọng của chỉ số chứng khoáng ACB sẽ tăng khoảng 0.17% với độ lệch chuẩn dự kiến sẽ là 1.48%. 3.2. Mô hình ARMAX–GARCH Sử dụng nguồn số liệu chỉ số chứng khoán ACB để phân tích mô hình, xem các ngày lễ có ảnh hưởng như thế nào đến biến động của thị trường chứng khoáng ta tiến hành kiểm định mô hình ARMAX – GARCH như sau: Ta tạo ra chuỗi dữ liệu suất sinh lợi với 2199 quan sát và biến giả dummy được hiểu là các ngày lễ lớn trong năm. Theo các dòng lệnh bên dưới để tính chuỗi tỷ suất sinh lợi và chuỗi biến giả dummy. > #Load dữ liệu > event<-read.csv(file.choose(),header=TRUE) > data<-read.csv(file.choose(),header=TRUE) > temp<-data[,c(2,6,15,16)] > close<-temp > for(i in 1:dim(close)[1]) close[i,]<-temp[length(data[,6])-i+1,] > #Tạo biến giả > dummy<-rep(0,dim(close)[1]) > for(i in 1:length(dummy)){ Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 60 + for(j in 1:dim(event)[1]){ + if(close[i,3]==event[j,1] & close[i,4]==event[j,2]) + dummy[i]<-1 + }} > dummy<-dummy[-1] > dt<-ts(close[,2],start=c(1,1),frequency=5) > ln.dt<-log(dt) > SSL<-diff(ln.dt) Sau khi thực hiện các dòng lệnh trên, chuỗi tỷ suất sinh lợi được gán trong đối tượng SSL và chuỗi biến giả được gán trong đối tượng dummy. Ta tìm một mô hình ARMAX cho chuỗi SSL với biến giả là dummy. Thực hiện các bước phân tích tương tự như mục 2.3.2 ta tìm được mô hình tối ưu là SARIMAX(0,0,1)x(0,0,1),s=5. > fit1<-arimax(SSL,order=c(0,0,1),sea=list(order=c(0,0,1),5),xreg=dummy) > res<-resid(fit1) > tsnew<-res^2 Thực hiện tương tự mục 3.1 ta thấy chuỗi có hiệu ứng ARCH và tìm được mô hình tối ưu là ARCH(8). > fit2<-garch(res,order=c(0,8),trace=0) Warning message: In sqrt(pred$e) : NaNs produced Dự báo Ngày tiếp theo của chuỗi là ngày 30/09/2015 không trùng vào ngày lễ nào trong năm nên giá trị biến giả sẽ là 0. > forecastGARCH(fitGARCH=fit2,fitARMA=fit1,trace=TRUE,r=7,newxreg=0) $ARCH a0 a1 a2 a3 a4 a5 5.57e-05 3.59e-01 2.64e-01 7.51e-02 6.82e-02 7.64e-02 a6 a7 a8 5.56e-02 1.70e-02 1.65e-01 $ARMA ma1 sma1 intercept xreg 0.151617 0.091441 -0.00031 0.001522 $forecast Point res res^2 SSL.forecast VAR.forecast 1 (439,2) 0.008095 8.21E-05 2 (439,3) -0.00038 0.000464 3 (439,4) -0.00347 6.84E-05 4 (439,5) 0.004703 6.11E-05 5 (440,1) 0.007817 2.21E-05 6 (440,2) 0.008268 1.20E-05 Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 61 7 (440,3) -0.02154 1.00E-07 8 (440,4) 0.009061 6.55E-05 9 (440,5) 0.0014507 0.00023 > sd<-sqrt( 0.0012795) > sd [1] 0.0357701 Theo kết quả mô hình ARCH(8), vào ngày 30/09/2015 không trùng vào ngày lễ nào trong năm nên ta có suất sinh lợi kỳ vọng của chỉ số chứng khoán ACB thay đổi không đáng kể (giảm khoảng 0.33%) với độ lệch chuẩn dự kiến sẽ là 3.58%. Phân tích chuỗi thời gian với sự hỗ trợ của package AnalyzeTS 62 TÀI LIỆU THAM KHẢO [1] Hồng Việt Minh, Luận văn tốt nghiệp Đại học: Phân tích số liệu thống kê với ngôn ngữ R. [2] Mai Thị Hồng Diễm, Luận văn tốt nghiệp Đại học: Phân tích chuỗi tài chính bằng mô hình chuỗi thời gian.