介紹

代碼、數據全部免費，都放在我的gitee倉庫里面https://gitee.com/yuanzhoulvpi/time_series，想要使用的可以直接到我這個倉庫下載。

本文是繼python與時間序列(開篇)詳細介紹了如何使用python計算時間序列模型

本文內容比較多【初略統計：所有文本超過1.9萬字，代碼有900行，49個圖】。

jupyter notebook目錄如下：

• 介紹時間序列

  • 導入時間序列數據

  • 處理(清洗)數據

  • 對數據進行可視化

  • 時間戳和時間段

  • 使用date_range

  • 滯后

  • 重采樣

  • 金融統計數據

  • 還手、回報、差分

  • 多個時間序列的比較

  • 窗口函數

  • OHLC曲線、燭箱圖

  • 自相關和偏自相關

• 時間序列的分解

  • 白噪聲

  • 隨機游走

• 使用python建模

  • AR模型

  • MA模型

  • ARMA、ARIMA、SARIMA模型

  • Var模型

  • 結構時間序列模型

  • 動態因子模型

介紹時間序列

導入需要的包

#?Importing?libraries import?os import?warnings warnings.filterwarnings('ignore') import?numpy?as?np import?pandas?as?pd import?matplotlib.pyplot?as?plt plt.style.use('fivethirtyeight') #?#?Above?is?a?special?style?template?for?matplotlib,?highly?useful?for?visualizing?time?series?data #?%matplotlib?inline #?from?pylab?import?rcParams #?from?plotly?import?tools #?import?plotly.plotly?as?py from?plotly.offline?import?init_notebook_mode,?iplot init_notebook_mode(connected=True) import?plotly.graph_objs?as?go import?plotly.figure_factory?as?ff import?statsmodels.api?as?sm from?numpy.random?import?normal,?seed from?scipy.stats?import?norm from?statsmodels.tsa.arima_model?import?ARMA from?statsmodels.tsa.stattools?import?adfuller from?statsmodels.graphics.tsaplots?import?plot_acf,?plot_pacf from?statsmodels.tsa.arima_process?import?ArmaProcess from?statsmodels.tsa.arima_model?import?ARIMA import?math import?seaborn?as?sns from?sklearn.metrics?import?mean_squared_error

導入時間序列數據

#?在使用pandas讀取文件的時候，可以使用index_col將數據的某一列直接設置為列索引； # parse_date這個參數可以將某一列數據自動解析為時間格式。google?=?pd.read_csv("../datasets/GOOGL_2006-01-01_to_2018-01-01.csv",?index_col='Date',?parse_dates=['Date']) google.head(5)#%%humidity?=?pd.read_csv("../datasets/humidity.csv",?index_col='datetime',?parse_dates=['datetime']) humidity.head(4)

處理(清洗)數據

上面的數據中, google數據是沒有缺失值的，但是humidity數據是有缺失值的，因此，需要對數據做一下處理。

#?觀察到humidity數據的第一行是空的，而且每一個時間基本上都非常接近的，因此我們將第一行給丟掉，然后使用向前填充數據humidity?=?humidity.iloc[1:] humidity?=?humidity.fillna(method='ffill') humidity.head(4)

對數據進行可視化

#?將時間序列轉換成特定的頻率，下面的.asfreq('M')中，就是將時間序列數據轉換成頻率為月的 fig,?ax?=?plt.subplots(figsize=(10,4),?dpi=300) humidity["Kansas?City"].asfreq('M').plot(ax=ax,?c='blue')?#?asfreq?method?is?used?to?convert?a?time?series?to?a?specified?frequency.?Here?it?is?monthly?frequency. plt.title('Humidity?in?Kansas?City?over?time(Monthly?frequency)') plt.show()bigcolor?=?['#1f77b4',?'#ff7f0e',?'#2ca02c',?'#d62728',?'#9467bd',?'#8c564b',?'#e377c2',?'#7f7f7f',?'#bcbd22',?'#17becf'] temp_colname?=?['Open','High',?'Low',?'Close',?'Volume'] fig,?ax?=?plt.subplots(nrows=len(temp_colname),?ncols=1,?figsize=(10,12),sharex=True,dpi=300)for?index,?colname_?in?enumerate(temp_colname):tempdata?=?google['2008':'2010'][colname_]ax[index].plot(tempdata.index,tempdata,?label=colname_,?c=bigcolor[index])ax[index].legend()#?google['2008':'2010'][colname_].plot(ax=ax[index],?label=colname_) fig.suptitle("Google?stock?attributes?from?2008?to?2010",y=0.92) plt.xticks(rotation=45)

時間戳和時間段

經常處理數據的可以遇到這樣的單詞：Timestamps，其實可以叫時間戳，本質上就是一個時間點。而Period,可以叫時間段，強調的是一段時間。

#%%#?創建一個時間點 timestamp?=?pd.Timestamp(2021,10,1,12) timestamp#%%#?創建一個時間段 period?=?pd.Period('2021-10-01') #?這個時間段，默認指的就是這一天 period#%%#?看看上面的時間點是不是在這個時間段里面 period.start_time?<?timestamp?<?period.end_time#%%#?將時間點轉換為時間段 new_period?=?timestamp.to_period(freq='H') new_period

使用date_range

date_range可以返回一個固定頻率的日期時間格式的索引，對已有的數據創建時間索引，然后處理數據。

#%%#?創建一個頻率以天的時間索引 dr1?=?pd.date_range(start='2021-10-01'?,end='2021-10-07') dr1#%%#?創建一個頻率為月的時間索引 dr2?=?pd.date_range(start='2020-10-01',?end='2021-11-01',?freq='M') dr2#%%#?創建一個以特定時間周期的時間索引（這里將一個時間段分成了5個時間點） dr3?=?pd.date_range(start='2021-10-01',?end='2021-11-14',?periods=5) dr3#%%#?如果只是同時設置start和period，會將時間向后衍生period天； dr4?=?pd.date_range(start='2021-10-01',periods=5) dr4#%%#?如果只是同時設置end和period，會將時間向前衍生period天； dr5?=?pd.date_range(end='2021-11-14',?periods=5) dr5

滯后

使用pandas的shift就行

#%%fig,?ax?=?plt.subplots(figsize=(10,4),?dpi=300) humidity['Vancouver'].asfreq('M').plot(ax=ax,?legend=True) shifted?=?humidity['Vancouver'].asfreq('M').shift(10).plot(ax=ax,?legend=True) shifted.legend(['Vancouver',?'Vancouver_lagged']) fig.suptitle(t="Vancouver_lagged_10?VS?Vancouver")

重采樣

上采樣：指的就是將數據從低頻數據上升到高頻數據，比如從月度數據上升到日數據。涉及到插入或者填充數據。

下采樣：指的就是將時間從高平數據下降到低頻數據，比如從日數據下降到月數據。涉及到的就是一些聚合操作。

#%%#?使用pressure數據來做這樣的操作 pressure?=?pd.read_csv('../datasets/pressure.csv',?index_col='datetime',?parse_dates=['datetime']) pressure.tail()#%%#?填充缺失數據，使用向前填充 pressure?=?pressure.iloc[1:] pressure?=?pressure.fillna(method='ffill') pressure.tail()#%%#?填充缺失數據，使用向后填充 pressure?=?pressure.fillna(method='bfill') pressure.head()#%%#?在采樣之前，查看pressure數據的大小 pressure.shape#%%#?開始降低采樣，3天一次 pressure?=?pressure.resample('3D').mean() pressure.head()#%%#?查看這個降采樣后的數據大小 pressure.shape#%%#?這個時候再進行上采樣 pressure?=?pressure.resample('D').pad() pressure.head()#%%#?上采樣之后，查看數據大小 pressure.shape

金融統計

查看換手率

fig,?ax?=?plt.subplots(figsize=(10,?4),?dpi=300) google['Change']?=?google.High.div(google.High.shift()) google['Change'].plot(ax=ax,?linewidth=1)

回報

#%%fig,?ax?=?plt.subplots(figsize=(10,?4),?dpi=300) google['Return']?=?google.Change.sub(1).mul(100) google['Return'].plot(ax=ax,?linewidth=1)#?另外一種方式計算回報fig,?ax?=?plt.subplots(figsize=(10,?4),?dpi=300) google.High.pct_change().mul(100).plot(ax=ax,?linewidth=1)

差分

fig,?ax?=?plt.subplots(figsize=(10,?4),?dpi=300) google.High.diff().plot(ax=ax,?linewidth=1)

比較兩個或多個時間序列數據

#?這里選擇微軟的數據和谷歌數據進行比較 microsoft?=?pd.read_csv('../datasets/MSFT_2006-01-01_to_2018-01-01.csv',?index_col='Date',?parse_dates=['Date'])#%%#?在還沒對數據進行標準化之前，看看原始數據是什么樣子的fig,?ax?=?plt.subplots(figsize=(10,?4),?dpi=300) google.High.plot(ax=ax) microsoft.High.plot(ax=ax) plt.legend(['Google','Microsoft'])#?對兩個數據進行標準化，然后再進行比較 #?將數據全部除以第一個時間點的數據，然后再乘以100 fig,?ax?=?plt.subplots(figsize=(10,?4),?dpi=300)normalized_google?=?google.High.div(google.High.iloc[0]).mul(100) normalized_microsoft?=?microsoft.High.div(microsoft.High.iloc[0]).mul(100) normalized_google.plot(ax=ax) normalized_microsoft.plot(ax=ax) plt.legend(['Google','Microsoft'])

基本上可以看出來，谷歌的增長要遠遠比微軟增長的要快一點

窗口函數

窗口函數可以識別子周期，可以計算子周期的一些性能數據

Rolling:窗口大小保持不變，窗口在數據上滑動，

Expanding: 窗口大小不斷拉長，計算的數據越來越多。

fig,?ax?=?plt.subplots(figsize=(10,?4),?dpi=300)#?Rolling?window?functions rolling_google?=?google.High.rolling('90D').mean() google.High.plot(ax=ax) rolling_google.plot(ax=ax) plt.legend(['High','Rolling?Mean'])

可以看出來，均值滾動的數據比原來的數據要更加平滑

fig,?ax?=?plt.subplots(figsize=(10,?4),?dpi=300)#?Expanding?window?functions microsoft_mean?=?microsoft.High.expanding().mean() microsoft_std?=?microsoft.High.expanding().std() microsoft.High.plot(ax=ax) microsoft_mean.plot(ax=ax) microsoft_std.plot(ax=ax) plt.legend(['High','Expanding?Mean','Expanding?Standard?Deviation'])

OHLC曲線

這個曲線說白了，就是顯示了股票數據的OPEN、HIGH、LOW、CLOSE的數據。如果一天的OPEN是大于CLOSE那么就是紅色的（因為當天跌了）；如果一天的OPEN是小于CLOSE就是綠色（因為當天漲了）

#?OHLC?chart?of?June?2008 trace?=?go.Ohlc(x=google['06-2008'].index,open=google['06-2008'].Open,high=google['06-2008'].High,low=google['06-2008'].Low,close=google['06-2008'].Close) #?data?=?[trace] #?iplot(data,?filename='simple_ohlc') fig?=?go.Figure() fig.add_trace(trace) fig.update_layout(template='plotly_white',title='sample?OHLC')#?OHLC?chart?of?2008 trace?=?go.Ohlc(x=google['2008'].index,open=google['2008'].Open,high=google['2008'].High,low=google['2008'].Low,close=google['2008'].Close) fig?=?go.Figure() fig.add_trace(trace) fig.update_layout(template='plotly_white',title='sample?OHLC')#%%#?OHLC?chart?of?2008 trace?=?go.Ohlc(x=google.index,open=google.Open,high=google.High,low=google.Low,close=google.Close) fig?=?go.Figure() fig.add_trace(trace) fig.update_layout(template='plotly_white',title='sample?OHLC')

燭箱圖

#?Candlestick?chart?of?march?2008 trace?=?go.Candlestick(x=google['03-2008'].index,open=google['03-2008'].Open,high=google['03-2008'].High,low=google['03-2008'].Low,close=google['03-2008'].Close) fig?=?go.Figure()fig.add_trace(trace) fig.update_layout(template='plotly_white',title='simple_candlestick')

自相關和偏自相關

#??humidity?of?San?Diego的自相關 fig,?ax?=?plt.subplots(figsize=(10,?4),?dpi=300) plot_acf(humidity["San?Diego"],lags=25,title="San?Diego",?ax=ax) plt.show()

由于所有滯后都接近 1 或至少大于置信區間，因此它們具有統計顯著性

#?humidity?of?San?Diego的偏自相關 fig,?ax?=?plt.subplots(figsize=(10,?4),?dpi=300) plot_pacf(humidity["San?Diego"],lags=25,ax=ax) plt.show()

在兩個滯后項以后，系數非常低了

時間序列的分解

任何時間序列可以可以被拆分為3個部分：

趨勢：趨勢是比較明顯的，比如極速的上升或者迅速下跌。

季節性：可以在數據中看到明顯的周期性，并且這個周期性和時間周期有關。這個周期可能是月，可能是季度，也可能是年。

誤差項。

但是不是說所有的時間序列都必須具備這3個部分。時間序列可能沒有明顯的趨勢、可能沒有明顯的周期性。或者兩個都沒有。

因此，可以將時間序列看成趨勢、周期性、誤差項的組合。

decomposed_google_volume?=?sm.tsa.seasonal_decompose(google["High"],freq=360)?#?The?frequncy?is?annual#?畫圖 fig,?ax?=?plt.subplots(ncols=1,?nrows=4,?figsize=(10,?12),?sharex=True) font?=?{'family':?'serif','color':?'darkred','weight':?'normal','size':?16,}def?plot_decompose(result,?ax,??title,?fontdict=font):ax[0].set_title(title,?fontdict=fontdict)result.observed.plot(ax=ax[0])ax[0].set_ylabel("Observed")result.trend.plot(ax=ax[1])ax[1].set_ylabel("Trend")result.seasonal.plot(ax=ax[2])ax[2].set_ylabel("Seasonal")result.resid.plot(ax=ax[3])ax[3].set_ylabel("Resid")plot_decompose(result=decomposed_google_volume,?ax=ax,?title="Google?High?decompose") plt.show()

可以看出來，谷歌的High有很強的趨勢性、季節性、而且上面的殘差還有一些模式信息沒有被提取完全

白噪聲

均值為固定；

方差固定；

滯后多少項，自相關都是0；

fig,?ax?=?plt.subplots(figsize=(10,?4),?dpi=300) white_noise?=?np.random.normal(loc=0,?scale=1,?size=1000) ax.plot(white_noise) plt.show()#?查看白噪聲的acffig,?ax?=?plt.subplots(ncols=1,nrows=1,figsize=(10,?4),?dpi=300) plot_acf(white_noise,lags=20,?ax=ax,?title="white_noise?acf") plt.show()

隨機游走

使用adf可以檢驗一個時間序列是否為隨機游走

#?Augmented?Dickey-Fuller?test?on?volume?of?google?and?microsoft?stocks adf?=?adfuller(microsoft["Volume"]) print("p-value?of?microsoft:?{}".format(float(adf[1]))) adf?=?adfuller(google["Volume"]) print("p-value?of?google:?{}".format(float(adf[1])))

上面兩個P值都是小于0.05，拒絕原假設，所以這兩個時間序列都不是隨機游走。

生成一個隨機游走數據

fig,?ax?=?plt.subplots(figsize=(10,?4),?dpi=300)random_walk?=?normal(loc=0,?scale=0.01,?size=1000) ax.plot(random_walk) plt.show()#?查看這個隨機游走模型的分布 fig,?ax?=?plt.subplots(figsize=(10,4),?dpi=300) sns.distplot(random_walk,?hist=True,?kde=True,bins=40,?color?=?'darkblue',hist_kws={'edgecolor':'black'},kde_kws={'linewidth':?4},?ax=ax) ax.set_title("Random?Walk")

使用統計工具箱

AR模型

自回歸模型是當前值可以被自身過去p階滯后的數據所解釋，p決定了需要幾個過去值來預測當前值。

AR(1)模型

如果上面的,這個時間蓄力就是隨機游走；如果，這個時間序列就是白噪聲；如果,這個時間序列就是平穩的。

AR(2)模型

AR(3)模型

模擬AR(1)模型

fig,?ax?=?plt.subplots(ncols=1,?nrows=4,?figsize=(10,?10),?dpi=300,?sharex=True) #?AR(1)?MA(1)?model:AR?parameter?=?+0.9ar1?=?np.array([1,?-0.9])?#?We?choose?-0.9?as?AR?parameter?is?+0.9 ma1?=?np.array([1]) AR1?=?ArmaProcess(ar1,?ma1) sim1?=?AR1.generate_sample(nsample=1000) ax[0].set_title('AR(1)?model:?AR?parameter?=?+0.9',fontsize=13) ax[0].plot(sim1)#?We?will?take?care?of?MA?model?later #?AR(1)?MA(1)?AR?parameter?=?-0.9 ar2?=?np.array([1,?0.9])?#?We?choose?+0.9?as?AR?parameter?is?-0.9 ma2?=?np.array([1]) AR2?=?ArmaProcess(ar2,?ma2) sim2?=?AR2.generate_sample(nsample=1000) ax[1].set_title('AR(1)?model:?AR?parameter?=?-0.9',fontsize=13) ax[1].plot(sim2)#?AR(2)?MA(1)?AR?parameter?=?0.9 plt.subplot(4,1,3) ar3?=?np.array([2,?-0.9])?#?We?choose?-0.9?as?AR?parameter?is?+0.9 ma3?=?np.array([1]) AR3?=?ArmaProcess(ar3,?ma3) sim3?=?AR3.generate_sample(nsample=1000) ax[2].set_title('AR(2)?model:?AR?parameter?=?+0.9',fontsize=13) ax[2].plot(sim3)#?AR(2)?MA(1)?AR?parameter?=?-0.9 plt.subplot(4,1,4) ar4?=?np.array([2,?0.9])?#?We?choose?+0.9?as?AR?parameter?is?-0.9 ma4?=?np.array([1]) AR4?=?ArmaProcess(ar4,?ma4) sim4?=?AR4.generate_sample(nsample=1000) ax[3].set_title('AR(2)?model:?AR?parameter?=?-0.9',fontsize=13) ax[3].plot(sim4)

對模擬數據做預測

model?=?ARMA(sim1,?order=(1,0)) result?=?model.fit() print(result.summary()) print("μ={}?,?={}".format(result.params[0],result.params[1]))

可以看出來，上面的?還是非常接近我們設置的模擬值的

#?將預測的數據畫出來 fig,?ax?=?plt.subplots(figsize=(10,?4),?dpi=300) #?Predicting?simulated?AR(1)?model result.plot_predict(start=900,?end=1010,ax=ax) plt.show()#?查看預測的rmse值 rmse?=?math.sqrt(mean_squared_error(sim1[900:1011],?result.predict(start=900,end=999))) print("The?root?mean?squared?error?is?{}.".format(rmse))fig,?ax?=?plt.subplots(figsize=(10,?3),?dpi=300)#?預測humidity?level?of?Montreal humid?=?ARMA(humidity["Montreal"].diff().iloc[1:].values,?order=(1,0)) res?=?humid.fit() res.plot_predict(start=1000,?end=1100,?ax=ax) plt.show()#?通過rmse查看預測效果 rmse?=?math.sqrt(mean_squared_error(humidity["Montreal"].diff().iloc[900:1000].values,?result.predict(start=900,end=999))) print("The?root?mean?squared?error?is?{}.".format(rmse))

上面的結果還不錯，我們再來繼續預測谷歌的數據

fig,?ax?=?plt.subplots(figsize=(10,?3),?dpi=300)#?Predicting?closing?prices?of?google humid?=?ARMA(google["Close"].diff().iloc[1:].values,?order=(1,0)) res?=?humid.fit() res.plot_predict(start=900,?end=1010,?ax=ax) plt.show()

MA模型

作為自回歸模型的替代，q階移動平均模型方程MA(q)中，左邊的是由右側的白噪聲線性組合得到的。

MA(1)模型

上面的公式表示，今天的值就是今天的噪聲加上昨天的一個噪聲乘以一個。

模擬MA(1)模型

fig,?ax?=?plt.subplots(figsize=(10,?3),?dpi=300) ar1?=?np.array([1]) ma1?=?np.array([1,?-0.5]) MA1?=?ArmaProcess(ar1,?ma1) sim1?=?MA1.generate_sample(nsample=1000) ax.plot(sim1) plt.show()

對模擬的MA數據做預測

model?=?ARMA(sim1,?order=(0,1)) result?=?model.fit() print(result.summary()) print("μ={}?,θ={}".format(result.params[0],result.params[1]))#?Forecasting?and?predicting?montreal?humidity model?=?ARMA(humidity["Montreal"].diff().iloc[1:].values,?order=(0,3)) result?=?model.fit() print(result.summary()) print("μ={}?,θ={}".format(result.params[0],result.params[1])) fig,?ax?=?plt.subplots(figsize=(10,?3),?dpi=300)result.plot_predict(start=1000,?end=1100,?ax=ax) plt.show()rmse?=?math.sqrt(mean_squared_error(humidity["Montreal"].diff().iloc[1000:1101].values,?result.predict(start=1000,end=1100))) print("The?root?mean?squared?error?is?{}.".format(rmse))

ARMA模型

我們選擇處理的是平穩時間序列中的自回歸、移動平均和混合自回歸移動平均(ARMA)模型的一般情況。

ARMA(p,q)模型公式如下：

使用ARMA模型做預測

#?Forecasting?and?predicting?microsoft?stocks?volume model?=?ARMA(microsoft["Volume"].diff().iloc[1:].values,?order=(3,3)) result?=?model.fit() print(result.summary()) print("μ={},??={},?θ={}".format(result.params[0],result.params[1],result.params[2])) fig,?ax?=?plt.subplots(figsize=(10,?3),?dpi=300)result.plot_predict(start=1000,?end=1100,ax=ax) plt.show()rmse?=?math.sqrt(mean_squared_error(microsoft["Volume"].diff().iloc[1000:1101].values,?result.predict(start=1000,end=1100))) print("The?root?mean?squared?error?is?{}.".format(rmse))

可以看出來，ARMA模型的效果比AR或者MA模型都要好的多

ARIMA模型

fig,?ax?=?plt.subplots(figsize=(10,?3),?dpi=300)#?Predicting?the?microsoft?stocks?volume model?=?ARIMA(microsoft["Volume"].diff().iloc[1:].values,?order=(2,1,0)) result?=?model.fit() print(result.summary()) result.plot_predict(start=700,?end=1000,?ax=ax) plt.show()rmse?=?math.sqrt(mean_squared_error(microsoft["Volume"].diff().iloc[700:1001].values,?result.predict(start=700,end=1000))) print("The?root?mean?squared?error?is?{}.".format(rmse))

Var模型

向量自回歸模型

#?Predicting?closing?price?of?Google?and?microsoft train_sample?=?pd.concat([google["Close"].diff().iloc[1:],microsoft["Close"].diff().iloc[1:]],axis=1) model?=?sm.tsa.VARMAX(train_sample,order=(2,1),trend='c') result?=?model.fit(maxiter=1000,disp=False) print(result.summary()) predicted_result?=?result.predict(start=0,?end=1000) #?fig,?ax?=?plt.subplots(figsize=(10,?10),?dpi=300) fig?=?plt.figure(figsize=(16,10),?dpi=300) result.plot_diagnostics(fig=fig) #?calculating?error rmse?=?math.sqrt(mean_squared_error(train_sample.iloc[1:1002].values,?predicted_result.values)) print("The?root?mean?squared?error?is?{}.".format(rmse))

SARIMA模型

時間序列季節自回歸求和滑動平均（SARIMA）

#?Predicting?closing?price?of?Google' train_sample?=?google["Close"].diff().iloc[1:].values model?=?sm.tsa.SARIMAX(train_sample,order=(4,0,4),trend='c') result?=?model.fit(maxiter=1000,disp=False) print(result.summary()) predicted_result?=?result.predict(start=0,?end=500) fig?=?plt.figure(figsize=(16,10),?dpi=300) result.plot_diagnostics(fig=fig) #?calculating?error rmse?=?math.sqrt(mean_squared_error(train_sample[1:502],?predicted_result)) print("The?root?mean?squared?error?is?{}.".format(rmse))

fig,?ax?=?plt.subplots(figsize=(10,?3),?dpi=300) ax.plot(train_sample[1:502],color='red') ax.plot(predicted_result,color='blue') ax.legend(['Actual','Predicted']) fig.suptitle('Google?Closing?prices') plt.show()

未觀察到的成分模型(又稱結構時間模型)

#?Predicting?closing?price?of?Google' train_sample?=?google["Close"].diff().iloc[1:].values model?=?sm.tsa.UnobservedComponents(train_sample,'local?level') result?=?model.fit(maxiter=1000,disp=False) print(result.summary()) predicted_result?=?result.predict(start=0,?end=500) fig?=?plt.figure(figsize=(16,10),?dpi=300)result.plot_diagnostics(fig=fig) #?calculating?error rmse?=?math.sqrt(mean_squared_error(train_sample[1:502],?predicted_result)) print("The?root?mean?squared?error?is?{}.".format(rmse))fig,?ax?=?plt.subplots(figsize=(10,?3),?dpi=300) ax.plot(train_sample[1:502],color='red') ax.plot(predicted_result,color='blue') ax.legend(['Actual','Predicted']) fig.suptitle('Google?Closing?prices') plt.show()

動態因子模型

#?Predicting?closing?price?of?Google?and?microsoft train_sample?=?pd.concat([google["Close"].diff().iloc[1:],microsoft["Close"].diff().iloc[1:]],axis=1) model?=?sm.tsa.DynamicFactor(train_sample,?k_factors=1,?factor_order=2) result?=?model.fit(maxiter=1000,disp=False) print(result.summary()) predicted_result?=?result.predict(start=0,?end=1000) fig?=?plt.figure(figsize=(16,10),?dpi=300) result.plot_diagnostics(fig=fig) #?calculating?error rmse?=?math.sqrt(mean_squared_error(train_sample.iloc[1:1002].values,?predicted_result.values)) print("The?root?mean?squared?error?is?{}.".format(rmse))

說明

本文文字參考的是：https://www.kaggle.com/thebrownviking20/everything-you-can-do-with-a-time-series#4.-Modelling-using-statstools

代碼和數據經過整理，存放在：https://gitee.com/yuanzhoulvpi/time_series

公眾號：AI蝸牛車保持謙遜、保持自律、保持進步個人微信備注：昵稱+學校/公司+方向 如果沒有備注不拉群！ 拉你進AI蝸牛車交流群

總結

以上是生活随笔為你收集整理的【时间序列】python与时间序列基本教程4(超过1.9万字 代码超过900行 包括49个图)...的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。