天氣資料是一個典型的時間序列資料，嘗試用arima來對天氣進行預測。

1、資料集：

採用的資料集為北京2011年1月1日到2020年3月23日的天氣情況。

採集工具：後裔採集器

網址：

http：//www。

tianqihoubao。com/lishi/

beijing。html

資料集中包含日期、天氣狀況、氣溫（最高氣溫、最低氣溫）、風力風向。

下文針對“最高溫度“進行分析。

2、利用r語言進行arima：

library（forecast）

library（tseries）

#1 load data

trainsize = 3370*0。9

testsize = 3370*0。1

data = read。csv（“/Users/jessica/PycharmProjects/machinelearning/dataanalysis/data/weather-beijing。csv”，nrows=trainsize）

testdata = read。csv（“/Users/jessica/PycharmProjects/machinelearning/dataanalysis/data/weather-beijing。csv”，skip=trainsize）

tempseries=ts（data［2］，frequency = 365）

plot。ts（tempseries）

tsdisplay（tempseries，xlab=“time”，ylab=“temperature”，lag。max = 40）

#2 對時間序列進行拆分，得到三個元件：seasonal，trend，random，可以用加法模型處理季節性

tempcomp <- decompose（tempseries）

plot（tempcomp）

#從圖中看出來，具有明顯的季節性因素。

#3 平穩性檢驗：

adf。test（tempseries）

# 透過Dickey-Fuller，可以看出序列平穩

# 取出季節性因素

tsdisplay（tempseries-tempcomp$seasonal，xlab=“time”，ylab=“temperature”，lag。max = 40）

# 非季節性p=1，q=9

# 對季節性進行差分

diff1 <- diff（tempseries，365）

plot（diff1）

tsdisplay（diff1，xlab=“time”，ylab=“diff1”，lag。max = 40）

# 非季節性p=1，q=9

# 確定引數方法1

# 根據acf和pacf圖形，arima（1，0，9）（0，1，0）［365］

# 確定引數方法2

auto。arima（tempseries，trace=T）

# 根據auto。arima，arima（1，0，3）（0，1，0）［365］

# arima模型 arima（x，order=，include。mean=，method=，transform。pars=，fixed=，seasonal=）

# -x：要進行模型擬合的序列名字。

# -order：指定模型階數。

# -include。mean：指定是否需要擬合常數項。

# -method：指定引數估計方法。

# -transform。pars：指定是否需要人為干預引數。

# -fixed：對疏係數模型指定疏係數的位置。

# -seasonal：指定季節模型的階數與季節週期，該選項的命令格式為：

# seasonal = list（order=c（P，D，Q），period = pi）

# （1）加法模型：P=0，Q=0

# （2）乘法模型：P，Q不全為零

# 方法一：使用引數（1，0，3）（0，1，0）［365］

fit1 <- Arima（tempseries， order=c（1，0，3），seasonal=list（order=c（0，1，0），period=365））

#結果：AIC=14581。06 AICc=14581。08 BIC=14610。5

fit1

tsdisplay（residuals（fit1）） #顯示殘差

#殘差白噪聲檢驗

for（i in 1：2） print（Box。test（fit1$residual，lag=6*i））

#殘差序列為白噪聲，故模型擬合正確

# forecast

f。p1<-forecast（fit1，h=testsize，level=c（99。5））

print（“mse is”）

sum（（f。p1［4］-testdata［2］）^2）/testsize

#mes is 19。51638

# 方法二：（1，0，9）（0，1，0）［365］

fit2 <- Arima（tempseries， order=c（1，0，9），seasonal=list（order=c（0，1，0），period=365））

#結果：AIC=14589。92 AICc=14590。02 BIC=14654。7

fit2

tsdisplay（residuals（fit2）） #顯示殘差

#殘差白噪聲檢驗

for（i in 1：2） print（Box。test（fit2$residual，lag=6*i））

#殘差序列為白噪聲，故模型擬合正確

# forecast

f。p2<-forecast（fit2，h=testsize，level=c（99。5））

print（“mse is”）

sum（（f。p2［4］-testdata［2］）^2）/testsize

#mse is 19。49468

3、結果

兩種方法中，效果較好的是（1，0，9）（0，1，0）［365］，此引數下mse is 19。49468。