Tuesday, September 23, 2014

Các hàm R phổ biến 9: Phân tích dãy số liệu theo thời gian (time series analysis)

Các hàm phổ biến cho sắp xếp dữ liệu theo thời gian. Các mô hình phân tích như AR, MA, và ARIMA cũng được mô tả qua. Các hàm này phần lớn đã có trong R, nhưng cũng có thể dùng package TTR cho các phân tích chuyên dụng. 


Phân tích 
Hàm R
Sắp xếp dữ liệu bằng hàm ts

births = scan("~/Documents/births.txt");
ts.births = ts(births, frequency=12, start=c(1946, 1))

Vẽ biểu đồ time series
plot.ts(ts.births)

Trung bình động (moving average) với bậc k
library(TTR);
SMA(ts.births, n=5)

Phân tích trend, season, và random: dữ liệu không có season

ma5 = SMA(ts.births, n=5)
plot.ts(ma5)

Phân tích trend, season, và random: dữ liệu có season

comp = decompose(ts.births);
comp;
plot(comp)

Điều chỉnh cho season
adj.births = ts.births-comp$seasonal;
plot(adj.births)

Tiên lượng bằng phương pháp simple exponential smoothing
library(forecast);
forecast = HoltWinters(ts.births, beta=F, gamma=F);
forecast;
plot(ts.births); plot(forecast)

Tiên lượng bằng phương pháp Holt exponential smoothing

forecast1 = HoltWinters(ts.births, gamma=F);
forecast1; plot(forecast1)

Tiên lượng bằng phương pháp Holt-Winters smoothing

forecast1 = HoltWinters(ts.births);
forecast1; plot(forecast1);

forecast2 = forecast.HoltWinters(forecast1, h=60)
plot.forecast(forecast2)

Mô hình ARIMA: sai phân (differencing), differences=k
skirts = scan("~/Documents/skirts.txt", skip=5);
ts.skirts = ts(skirts, start=c(1866));
# tính sai phân
dif.skirts = diff(ts.skirts, differences=1);
plot.ts(dif.skirts)

Autocorrelation và partial autocorrelation

acf(dif.kings, lag.max=20, plot=F);  
acf(dif.kings, lag.max=20);
pacf(dif.kings, lag.max=20, plot=F); 
pacf(dif.kings, lag.max=20)


Mô hình MA
kings = scan("~/Documents/kings.txt")
library(forecast);
auto.arima(kings)

Mô hình RA
vc = scan("~/Documents/vulcano.txt")
ts.vc = ts(vc, start=c(1500));
plot.ts(ts.vc);
auto.arima(vc);
auto.arima(vc, ic="bic")

Mô hình ARIMA(p, d, q)
arima1 = arima(ts.kings, order=c(0, 1, 1));
arima1

Tiên lượng bằng mô hình ARIMA(p, d, q)
library(forecast);
forecast1 = forecast.Arima(arima1, h=5);
forecast1;
plot.forecast(forecast1)




No comments:

Post a Comment