十大ppt模板免費下載網(wǎng)站,wordpress遠程后臺設(shè)置,制作個人網(wǎng)頁詳細教程,最新開的手游傳奇網(wǎng)站時間序列分析#xff1a;從隨機過程理論到 SARIMA 模型構(gòu)建 一、 時間序列的數(shù)學(xué)定義與性質(zhì) 1.1 隨機序列與觀察值 在數(shù)學(xué)上#xff0c;時間序列不僅僅是一組數(shù)字#xff0c;它是一個隨機過程#xff08;Stochastic Process#xff09;。 設(shè) TTT 為一個時間索引集合#…時間序列分析從隨機過程理論到 SARIMA 模型構(gòu)建一、 時間序列的數(shù)學(xué)定義與性質(zhì)1.1 隨機序列與觀察值在數(shù)學(xué)上時間序列不僅僅是一組數(shù)字它是一個隨機過程Stochastic Process。設(shè)TTT為一個時間索引集合通常為整數(shù)集Zmathbb{Z}Z我們稱隨機變量族{Xt:t∈T}{X_t : t in T}{Xt?:t∈T}為一個時間序列。我們在現(xiàn)實中記錄到的具體數(shù)據(jù)如某年的GDP是這個隨機過程的一個樣本實現(xiàn)Realization記為{xt:t1,2,…,n}{x_t : t1, 2, dots, n}{xt?:t1,2,…,n}。importpandasaspdimportnumpyasnpimportmatplotlib.pyplotaspltimportseabornassns# 設(shè)置繪圖風(fēng)格sns.set(stylewhitegrid)plt.rcParams[font.sans-serif][SimHei]# 用來正常顯示中文標簽plt.rcParams[axes.unicode_minus]False# 用來正常顯示負號# 讀取數(shù)據(jù)# 假設(shè) Excel 文件名為 export_datas.xlsx時間列為 date數(shù)據(jù)列為 exportdfpd.read_excel(export_datas.xlsx)# 創(chuàng)建時間索引 (2001年1月 - 2016年8月)df[date]pd.date_range(start2001-01-01,periodslen(df),freqM)df.set_index(date,inplaceTrue)# 繪制原始時序圖plt.figure(figsize(12,6))plt.plot(df.index,df[export],label出口額(億美元),colornavy)plt.title(我國出口額月度數(shù)據(jù) (2001-2016),fontsize14)plt.xlabel(年份)plt.ylabel(出口額)plt.legend()plt.show()1.2 平穩(wěn)性 (Stationarity)平穩(wěn)性是時間序列建模的基石。如果一個隨機過程的統(tǒng)計規(guī)律不隨時間推移而改變則稱其為平穩(wěn)的。1. 嚴平穩(wěn) (Strict Stationarity)若對任意的t1,t2,…,tk∈Tt_1, t_2, dots, t_k in Tt1?,t2?,…,tk?∈T和任意整數(shù)τ auτ隨機變量組(Xt1,…,Xtk)(X_{t_1}, dots, X_{t_k})(Xt1??,…,Xtk??)與(Xt1τ,…,Xtkτ)(X_{t_1 au}, dots, X_{t_k au})(Xt1?τ?,…,Xtk?τ?)具有相同的聯(lián)合概率分布即Ft1,…,tk(x1,…,xk)Ft1τ,…,tkτ(x1,…,xk) F_{t_1, dots, t_k}(x_1, dots, x_k) F_{t_1 au, dots, t_k au}(x_1, dots, x_k)Ft1?,…,tk??(x1?,…,xk?)Ft1?τ,…,tk?τ?(x1?,…,xk?)則稱該序列為嚴平穩(wěn)序列。條件極其苛刻實際應(yīng)用中很少滿足。2. 寬平穩(wěn) (Weak Stationarity)在實際建模中我們通常使用寬平穩(wěn)或稱二階平穩(wěn)。一個時間序列{Xt}{X_t}{Xt?}稱為寬平穩(wěn)的若滿足以下三個條件均值恒定對于任意ttt期望存在且為常數(shù)。E(Xt)μ,?t∈T E(X_t) mu, quad forall t in TE(Xt?)μ,?t∈T二階矩有限E(Xt2)∞,?t∈T E(X_t^2) infty, quad forall t in TE(Xt2?)∞,?t∈T自協(xié)方差僅依賴于時間間隔對于任意t,st, st,s協(xié)方差僅與時間差ks?tk s-tks?t有關(guān)。Cov(Xt,Xtk)E[(Xt?μ)(Xtk?μ)]γk Cov(X_t, X_{tk}) E[(X_t - mu)(X_{tk} - mu)] gamma_kCov(Xt?,Xtk?)E[(Xt??μ)(Xtk??μ)]γk?二、 滯后算子與差分為了簡化模型的數(shù)學(xué)表達引入**滯后算子Lag Operator,LLL和差分算子Difference Operator,ΔDeltaΔ**是必要的。2.1 滯后算子LLL滯后算子作用于時間序列使其時間標號后退一期LXtXt?1 L X_t X_{t-1}LXt?Xt?1?推廣到kkk階LkXtXt?k L^k X_t X_{t-k}LkXt?Xt?k?常數(shù)ccc在滯后算子作用下不變LccL c cLcc。2.2 差分運算差分是處理非平穩(wěn)序列去除趨勢的主要手段。一階差分ΔXtXt?Xt?1(1?L)Xt Delta X_t X_t - X_{t-1} (1 - L)X_tΔXt?Xt??Xt?1?(1?L)Xt?kkk步差分季節(jié)差分ΔkXtXt?Xt?k(1?Lk)Xt Delta_k X_t X_t - X_{t-k} (1 - L^k)X_tΔk?Xt?Xt??Xt?k?(1?Lk)Xt?ddd階差分連續(xù)進行ddd次一階差分ΔdXt(1?L)dXt Delta^d X_t (1 - L)^d X_tΔdXt?(1?L)dXt?利用二項式定理展開其通項公式為ΔdXt∑j0d(?1)j(dj)Xt?j Delta^d X_t sum_{j0}^d (-1)^j inom0ghprlcq{j} X_{t-j}ΔdXt?j0∑d?(?1)j(jd?)Xt?j?# 1. 對數(shù)變換 (穩(wěn)定方差)df[log_export]np.log(df[export])# 2. 一階差分 (消除趨勢)# 注意差分后會產(chǎn)生 NaN 值需要去除df[diff_log_export]df[log_export].diff().dropna()# 繪制處理后的數(shù)據(jù)fig,axplt.subplots(2,1,figsize(12,10))# 對數(shù)數(shù)據(jù)圖ax[0].plot(df.index,df[log_export],colorgreen)ax[0].set_title(對數(shù)變換后的序列)# 差分數(shù)據(jù)圖ax[1].plot(df.index,df[diff_log_export],colororange)ax[1].set_title(對數(shù)一階差分后的序列)plt.tight_layout()plt.show()平穩(wěn)性檢驗fromstatsmodels.tsa.stattoolsimportadfullerdefadf_test(series): 封裝 ADF 檢驗函數(shù)輸出清晰的結(jié)果 print(ADF 檢驗結(jié)果:)resultadfuller(series.dropna())print(f統(tǒng)計量 (ADF Statistic):{result[0]:.4f})print(fp值 (p-value):{result[1]:.4f})print(臨界值 (Critical Values):)forkey,valueinresult[4].items():print(f{key}:{value:.4f})ifresult[1]0.05:print(結(jié)論: 拒絕原假設(shè)序列是平穩(wěn)的)else:print(結(jié)論: 接受原假設(shè)序列是非平穩(wěn)的)# 對差分后的序列進行檢驗adf_test(df[diff_log_export])三、 隨機過程ARIMA 模型是由以下幾個基礎(chǔ)過程構(gòu)建而成的。3.1 白噪聲過程 (White Noise)白噪聲是時間序列建模中的基本“原子”通常記為{εt}{varepsilon_t}{εt?}。若滿足以下性質(zhì)則稱為純隨機過程或白噪聲E(εt)0 E(varepsilon_t) 0E(εt?)0Var(εt)σ2∞ Var(varepsilon_t) sigma^2 inftyVar(εt?)σ2∞Cov(εt,εtk)0,?k≠0 Cov(varepsilon_t, varepsilon_{tk}) 0, quad forall k eq 0Cov(εt?,εtk?)0,?k03.2 自回歸過程 AR§ppp階自回歸模型 (AutoRegressive model of orderppp)描述了當(dāng)前值與過去ppp個歷史值之間的線性關(guān)系。數(shù)學(xué)表達式Xtc?1Xt?1?2Xt?2??pXt?pεt X_t c phi_1 X_{t-1} phi_2 X_{t-2} dots phi_p X_{t-p} varepsilon_tXt?c?1?Xt?1??2?Xt?2???p?Xt?p?εt?其中ccc為常數(shù)項?iphi_i?i?為自回歸系數(shù)εtvarepsilon_tεt?為白噪聲。引入滯后算子多項式Φ(L)Phi(L)Φ(L)令Φ(L)1??1L??2L2????pLpPhi(L) 1 - phi_1 L - phi_2 L^2 - dots - phi_p L^pΦ(L)1??1?L??2?L2????p?Lp則模型可簡寫為Φ(L)Xtcεt Phi(L)X_t c varepsilon_tΦ(L)Xt?cεt?平穩(wěn)性條件AR§ 序列平穩(wěn)的充要條件是特征方程Φ(z)0Phi(z) 0Φ(z)0的所有根zzz都在單位圓外即∣z∣1,for all z such that 1?∑i1p?izi0 |z| 1, quad ext{for all } z ext{ such that } 1 - sum_{i1}^p phi_i z^i 0∣z∣1,for allzsuch that1?i1∑p??i?zi03.3 移動平均過程 MA(q)qqq階移動平均模型 (Moving Average model of orderqqq)描述了當(dāng)前值與過去qqq個隨機擾動項誤差之間的線性關(guān)系。數(shù)學(xué)表達式Xtμεt?θ1εt?1?θ2εt?2???θqεt?q X_t mu varepsilon_t - heta_1 varepsilon_{t-1} - heta_2 varepsilon_{t-2} - dots - heta_q varepsilon_{t-q}Xt?μεt??θ1?εt?1??θ2?εt?2????θq?εt?q?(注根據(jù) Box-Jenkins 習(xí)慣MA 參數(shù)前常取負號也有教材取正號本質(zhì)只需調(diào)整參數(shù)符號即可)引入滯后算子多項式Θ(L)Theta(L)Θ(L)令Θ(L)1?θ1L?θ2L2???θqLqTheta(L) 1 - heta_1 L - heta_2 L^2 - dots - heta_q L^qΘ(L)1?θ1?L?θ2?L2???θq?Lq則模型可簡寫為XtμΘ(L)εt X_t mu Theta(L)varepsilon_tXt?μΘ(L)εt?性質(zhì)MA(q) 模型總是平穩(wěn)的。但為了保證模型的可逆性 (Invertibility)即可以用過去的觀察值表示當(dāng)前的誤差要求方程Θ(z)0Theta(z) 0Θ(z)0的所有根都在單位圓外。3.4 自回歸移動平均模型 ARMA(p, q)將 AR 和 MA 結(jié)合得到 ARMA 模型Xtc∑i1p?iXt?iεt?∑j1qθjεt?j X_t c sum_{i1}^p phi_i X_{t-i} varepsilon_t - sum_{j1}^q heta_j varepsilon_{t-j}Xt?ci1∑p??i?Xt?i?εt??j1∑q?θj?εt?j?算子形式Φ(L)XtcΘ(L)εt Phi(L)X_t c Theta(L)varepsilon_tΦ(L)Xt?cΘ(L)εt?fromstatsmodels.graphics.tsaplotsimportplot_acf,plot_pacf# 繪制 ACF 和 PACF 圖fig,axplt.subplots(2,1,figsize(12,8))# 自相關(guān)圖 (ACF)plot_acf(df[diff_log_export].dropna(),lags40,axax[0])ax[0].set_title(自相關(guān)函數(shù) (ACF))# 偏自相關(guān)圖 (PACF)plot_pacf(df[diff_log_export].dropna(),lags40,axax[1])ax[1].set_title(偏自相關(guān)函數(shù) (PACF))plt.tight_layout()plt.show()# 提示也可以使用 pmdarima 庫進行自動定階# import pmdarima as pm# model pm.auto_arima(df[log_export], seasonalTrue, m12)四、 ARIMA 與 SARIMA非平穩(wěn)與季節(jié)性建模4.1 ARIMA(p, d, q) 模型現(xiàn)實數(shù)據(jù)往往是非平穩(wěn)的含有趨勢。如果一個時間序列{Yt}{Y_t}{Yt?}經(jīng)過ddd次差分后變?yōu)槠椒€(wěn)的 ARMA(p, q) 過程則稱{Yt}{Y_t}{Yt?}服從 ARIMA(p, d, q) 模型。令XtΔdYt(1?L)dYtX_t Delta^d Y_t (1-L)^d Y_tXt?ΔdYt?(1?L)dYt?且XtX_tXt?滿足Φ(L)XtΘ(L)εtPhi(L)X_t Theta(L)varepsilon_tΦ(L)Xt?Θ(L)εt?。ARIMA 模型的完整算子公式Φ(L)(1?L)dYtcΘ(L)εt Phi(L) (1 - L)^d Y_t c Theta(L) varepsilon_tΦ(L)(1?L)dYt?cΘ(L)εt?展開即為(1?∑i1p?iLi)(1?L)dYtc(1?∑j1qθjLj)εt left(1 - sum_{i1}^p phi_i L^i ight) (1 - L)^d Y_t c left(1 - sum_{j1}^q heta_j L^j ight) varepsilon_t(1?i1∑p??i?Li)(1?L)dYt?c(1?j1∑q?θj?Lj)εt?4.2 SARIMA 模型季節(jié)性的引入對于具有周期性如月度周期S12S12S12季度周期S4S4S4的數(shù)據(jù)普通的 ARIMA 無法捕捉季節(jié)效應(yīng)。我們需要引入季節(jié)性自回歸移動平均模型 (SARIMA)。SARIMA 模型的結(jié)構(gòu)由非季節(jié)部分(p,d,q)(p, d, q)(p,d,q)和季節(jié)部分(P,D,Q)S(P, D, Q)_S(P,D,Q)S?組成記為ARIMA(p,d,q)×(P,D,Q)S ARIMA(p, d, q) imes (P, D, Q)_SARIMA(p,d,q)×(P,D,Q)S?1. 季節(jié)算子定義季節(jié)性 AR 多項式(階數(shù)PPP)ΦS(LS)1?Φ1LS?Φ2L2S???ΦPLPS Phi_S(L^S) 1 - Phi_1 L^S - Phi_2 L^{2S} - dots - Phi_P L^{PS}ΦS?(LS)1?Φ1?LS?Φ2?L2S???ΦP?LPS季節(jié)性 MA 多項式(階數(shù)QQQ)ΘS(LS)1?Θ1LS?Θ2L2S???ΘQLQS Theta_S(L^S) 1 - Theta_1 L^S - Theta_2 L^{2S} - dots - Theta_Q L^{QS}ΘS?(LS)1?Θ1?LS?Θ2?L2S???ΘQ?LQS季節(jié)性差分算子(階數(shù)DDD)ΔSD(1?LS)D Delta_S^D (1 - L^S)^DΔSD?(1?LS)D2. SARIMA 通用數(shù)學(xué)表達式Φ(L)ΦS(LS)ΔdΔSDYtΘ(L)ΘS(LS)εt Phi(L) Phi_S(L^S) Delta^d Delta_S^D Y_t Theta(L) Theta_S(L^S) varepsilon_tΦ(L)ΦS?(LS)ΔdΔSD?Yt?Θ(L)ΘS?(LS)εt?3. 公式的詳細結(jié)構(gòu)解析為了完整展示該模型我們將上述各項代入并展開(1?∑i1p?iLi)?非季節(jié) AR(1?∑j1PΦjLjS)?季節(jié) AR(1?L)d?差分(1?LS)D?季節(jié)差分Yt(1?∑k1qθkLk)?非季節(jié) MA(1?∑m1QΘmLmS)?季節(jié) MAεt underbrace{left(1 - sum_{i1}^p phi_i L^i ight)}_{ ext{非季節(jié) AR}} underbrace{left(1 - sum_{j1}^P Phi_j L^{jS} ight)}_{ ext{季節(jié) AR}} underbrace{(1-L)^d}_{ ext{差分}} underbrace{(1-L^S)^D}_{ ext{季節(jié)差分}} Y_t underbrace{left(1 - sum_{k1}^q heta_k L^k ight)}_{ ext{非季節(jié) MA}} underbrace{left(1 - sum_{m1}^Q Theta_m L^{mS} ight)}_{ ext{季節(jié) MA}} varepsilon_t非季節(jié)AR(1?i1∑p??i?Li)??季節(jié)AR(1?j1∑P?Φj?LjS)??差分(1?L)d??季節(jié)差分(1?LS)D??Yt?非季節(jié)MA(1?k1∑q?θk?Lk)??季節(jié)MA(1?m1∑Q?Θm?LmS)??εt?這是時間序列分析中最為復(fù)雜也最為通用的線性模型形式。importstatsmodels.apiassm# 定義 SARIMA 模型# order(p, d, q) - (2, 1, 1)# seasonal_order(P, D, Q, s) - (2, 1, 1, 12)modelsm.tsa.statespace.SARIMAX(df[log_export],order(2,1,1),seasonal_order(2,1,1,12),enforce_stationarityFalse,enforce_invertibilityFalse)# 擬合模型resultsmodel.fit()# 打印模型詳細統(tǒng)計摘要print(results.summary())五、 參數(shù)估計與模型選擇5.1 極大似然估計 (MLE)在確定了模型階數(shù)(p,d,q)(p, d, q)(p,d,q)和(P,D,Q)(P, D, Q)(P,D,Q)后需要估計參數(shù)向量β(?,θ,Φ,Θ,σ2)eta (phi, heta, Phi, Theta, sigma^2)β(?,θ,Φ,Θ,σ2)。通常假設(shè)殘差服從正態(tài)分布εt～N(0,σ2)varepsilon_t sim N(0, sigma^2)εt?～N(0,σ2)使用對數(shù)似然函數(shù)ln?L(β∣Y)≈?n2ln?(2π)?n2ln?(σ2)?12σ2∑t1nεt2(β) ln L(eta | Y) approx -frac{n}{2}ln(2pi) - frac{n}{2}ln(sigma^2) - frac{1}{2sigma^2} sum_{t1}^n varepsilon_t^2(eta)lnL(β∣Y)≈?2n?ln(2π)?2n?ln(σ2)?2σ21?t1∑n?εt2?(β)通過數(shù)值優(yōu)化算法如 BFGS 或 Newton-Raphson最大化該函數(shù)以求解參數(shù)。5.2 信息準則 (Information Criteria)為了在模型擬合優(yōu)度與模型復(fù)雜度之間取得平衡避免過擬合我們使用信息準則進行模型選擇。1. 赤池信息準則 (AIC)AIC?2ln?(L^)2k AIC -2 ln(hat{L}) 2kAIC?2ln(L^)2k其中L^hat{L}L^是似然函數(shù)的最大值kkk是模型中被估計參數(shù)的總數(shù)kpqPQ1k p q P Q 1kpqPQ1。2. 貝葉斯信息準則 (BIC/SBC)BIC?2ln?(L^)kln?(n) BIC -2 ln(hat{L}) k ln(n)BIC?2ln(L^)kln(n)其中nnn為樣本容量。BIC 對參數(shù)個數(shù)的懲罰力度比 AIC 更大當(dāng)n≥8n geq 8n≥8時傾向于選擇更簡單的模型。判別規(guī)則AIC 或 BIC 值越小模型越優(yōu)。六、 殘差檢驗建模的最后一步是驗證模型假設(shè)即檢驗殘差序列{ε^t}{hat{varepsilon}_t}{ε^t?}是否為白噪聲。如果殘差中仍包含自相關(guān)性說明模型未能提取全部信息。Ljung-Box Q 檢驗原假設(shè)H0H_0H0?殘差序列不存在自相關(guān)即ρ1ρ2?ρm0 ho_1 ho_2 dots ho_m 0ρ1?ρ2??ρm?0。檢驗統(tǒng)計量QQQ定義為Qn(n2)∑k1mρ^k2n?k Q n(n2) sum_{k1}^m frac{hat{ ho}_k^2}{n-k}Qn(n2)k1∑m?n?kρ^?k2??其中nnn是樣本量。ρ^khat{ ho}_kρ^?k?是殘差的kkk階樣本自相關(guān)系數(shù)ρ^k∑tk1nε^tε^t?k∑t1nε^t2 hat{ ho}_k frac{sum_{tk1}^n hat{varepsilon}_t hat{varepsilon}_{t-k}}{sum_{t1}^n hat{varepsilon}_t^2}ρ^?k?∑t1n?ε^t2?∑tk1n?ε^t?ε^t?k??mmm是檢驗的最大滯后階數(shù)。統(tǒng)計量QQQ近似服從自由度為m?(pqPQ)m - (pqPQ)m?(pqPQ)的卡方分布Q～χ2(m?params) Q sim chi^2(m - ext{params})Q～χ2(m?params)若計算出的ppp值小于顯著性水平αalphaα如 0.05則拒絕原假設(shè)認為模型不合格反之則認為模型有效。fromstatsmodels.stats.diagnosticimportacorr_ljungbox# 1. 繪制殘差診斷圖 (包括直方圖、QQ圖、相關(guān)圖)results.plot_diagnostics(figsize(15,12))plt.show()# 2. Ljung-Box 白噪聲檢驗# 檢驗殘差是否存在自相關(guān)選取滯后 6 期和 12 期lb_testacorr_ljungbox(results.resid,lags[6,12],return_dfTrue)print(Ljung-Box 檢驗結(jié)果:)print(lb_test)# 如果 lb_pvalue 0.05則說明殘差為白噪聲模型通過檢驗?zāi)Ｐ皖A(yù)測# 預(yù)測未來 10 個月forecast_steps10pred_ucresults.get_forecast(stepsforecast_steps)# 獲取預(yù)測的均值和置信區(qū)間pred_mean_logpred_uc.predicted_mean pred_ci_logpred_uc.conf_int()# *** 關(guān)鍵步驟逆變換 (從對數(shù)還原為原始值) ***pred_meannp.exp(pred_mean_log)pred_cinp.exp(pred_ci_log)# 創(chuàng)建未來時間索引future_datespd.date_range(startdf.index[-1]pd.DateOffset(months1),periodsforecast_steps,freqM)pred_mean.indexfuture_dates pred_ci.indexfuture_dates# 打印預(yù)測結(jié)果forecast_dfpd.DataFrame({預(yù)測值:pred_mean,下限 (95%):pred_ci.iloc[:,0],上限 (95%):pred_ci.iloc[:,1]})print(forecast_df)# 繪圖展示plt.figure(figsize(14,7))# 繪制歷史數(shù)據(jù)plt.plot(df.index,df[export],label歷史觀察值,colorblack)# 繪制預(yù)測值plt.plot(pred_mean.index,pred_mean,label未來預(yù)測值,colorred,markero)# 繪制置信區(qū)間plt.fill_between(pred_ci.index,pred_ci.iloc[:,0],pred_ci.iloc[:,1],colorpink,alpha0.3,label95% 置信區(qū)間)plt.title(SARIMA 模型未來 10 個月出口額預(yù)測,fontsize16)plt.legend(locupper left)plt.grid(True)plt.show()七、 案例數(shù)學(xué)模型假設(shè)我們最終確定的模型為ARIMA(2,1,1)×(2,1,1)12ARIMA(2,1,1) imes (2,1,1)_{12}ARIMA(2,1,1)×(2,1,1)12?且對原始數(shù)據(jù)YtY_tYt?進行了對數(shù)變換ytln?(Yt)y_t ln(Y_t)yt?ln(Yt?)。這意味著我們需要對yty_tyt?進行一次普通差分(1?L)(1-L)(1?L)和一次周期為 12 的季節(jié)差分(1?L12)(1-L^{12})(1?L12)。該模型的具體數(shù)學(xué)方程展開如下(1??1L??2L2)(1?Φ1L12?Φ2L24)(1?L)(1?L12)ln?(Yt)(1?θ1L)(1?Θ1L12)εt (1 - phi_1 L - phi_2 L^2)(1 - Phi_1 L^{12} - Phi_2 L^{24}) (1-L)(1-L^{12}) ln(Y_t) (1 - heta_1 L)(1 - Theta_1 L^{12}) varepsilon_t(1??1?L??2?L2)(1?Φ1?L12?Φ2?L24)(1?L)(1?L12)ln(Yt?)(1?θ1?L)(1?Θ1?L12)εt?模型含義解析左邊ln?(Yt)ln(Y_t)ln(Yt?)對數(shù)變換解決異方差。(1?L)(1?L12)(1-L)(1-L^{12})(1?L)(1?L12)雙重差分消除長期趨勢和季節(jié)性趨勢。(1??1L??2L2)(1 - phi_1 L - phi_2 L^2)(1??1?L??2?L2)當(dāng)前的差分值與過去 1 個月和 2 個月的差分值相關(guān)。(1?Φ1L12?Φ2L24)(1 - Phi_1 L^{12} - Phi_2 L^{24})(1?Φ1?L12?Φ2?L24)當(dāng)前的季節(jié)性變化與去年的同月、前年的同月的變化相關(guān)。右邊(1?θ1L)(1 - heta_1 L)(1?θ1?L)包含 1 階非季節(jié)性移動平均成分短期震蕩。(1?Θ1L12)(1 - Theta_1 L^{12})(1?Θ1?L12)包含 1 階季節(jié)性移動平均成分季節(jié)性震蕩。εtvarepsilon_tεt?剩余的純隨機波動白噪聲。通過求解上述方程的期望即可得到對未來時刻YnhY_{nh}Ynh?的最優(yōu)預(yù)測值。八、完整代碼# -*- coding: utf-8 -*- 時間序列分析完整實戰(zhàn)代碼 (SARIMA模型) 對應(yīng)博客章節(jié)七、 Python 實戰(zhàn)我國出口額月度數(shù)據(jù)預(yù)測 importpandasaspdimportnumpyasnpimportmatplotlib.pyplotaspltimportseabornassnsimportstatsmodels.apiassmfromstatsmodels.tsa.stattoolsimportadfullerfromstatsmodels.graphics.tsaplotsimportplot_acf,plot_pacffromstatsmodels.stats.diagnosticimportacorr_ljungboximportwarnings# 忽略一些不必要的警告warnings.filterwarnings(ignore)# # 0. 全局設(shè)置# # 設(shè)置繪圖風(fēng)格sns.set(stylewhitegrid)# 設(shè)置字體以支持中文顯示# Windows用戶通常使用 SimHeiMac用戶使用 Arial Unicode MSplt.rcParams[font.sans-serif][SimHei]plt.rcParams[axes.unicode_minus]False# 開關(guān)是否使用模擬數(shù)據(jù)# 如果你有真實的 export_datas.xlsx 文件請將此處改為 FalseUSE_MOCK_DATATrue# # Step 1: 讀取數(shù)據(jù)與初步觀察# print( Step 1: 數(shù)據(jù)準備中...)ifUSE_MOCK_DATA:# --- 生成模擬數(shù)據(jù) (為了演示代碼可運行) ---datespd.date_range(start2001-01-01,periods188,freqM)# 模擬趨勢 季節(jié)性 隨機噪聲np.random.seed(42)trendnp.linspace(20,200,188)seasonality10*np.sin(np.linspace(0,3.14*16,188))noisenp.random.normal(0,5,188)data_valuestrendseasonalitynoise100# 保證為正數(shù)dfpd.DataFrame({date:dates,export:data_values})df.set_index(date,inplaceTrue)print(提示正在使用模擬數(shù)據(jù)運行...)else:# --- 讀取真實數(shù)據(jù) ---try:dfpd.read_excel(export_datas.xlsx)# 假設(shè)Excel中有 date 和 export 兩列ifdatenotindf.columns:# 如果沒有日期列生成一個df[date]pd.date_range(start2001-01-01,periodslen(df),freqM)df.set_index(date,inplaceTrue)exceptFileNotFoundError:print(錯誤未找到 export_datas.xlsx 文件。請將 USE_MOCK_DATA 設(shè)置為 True 以測試代碼。)exit()# 繪制原始時序圖plt.figure(figsize(12,6))plt.plot(df.index,df[export],label出口額(億美元),colornavy)plt.title(我國出口額月度數(shù)據(jù) (2001-2016),fontsize14)plt.xlabel(年份)plt.ylabel(出口額)plt.legend()plt.show()# # Step 2: 數(shù)據(jù)平穩(wěn)化處理# print( Step 2: 進行對數(shù)變換與差分處理...)# 1. 對數(shù)變換 (穩(wěn)定方差)df[log_export]np.log(df[export])# 2. 一階差分 (消除趨勢)df[diff_log_export]df[log_export].diff()# 繪制處理后的數(shù)據(jù)對比圖fig,axplt.subplots(2,1,figsize(12,10))ax[0].plot(df.index,df[log_export],colorgreen)ax[0].set_title(對數(shù)變換后的序列)ax[0].grid(True)ax[1].plot(df.index,df[diff_log_export],colororange)ax[1].set_title(對數(shù)一階差分后的序列)ax[1].grid(True)plt.tight_layout()plt.show()# # Step 3: 平穩(wěn)性檢驗 (ADF Test)# print( Step 3: 進行 ADF 平穩(wěn)性檢驗...)defadf_test(series):ADF 檢驗封裝函數(shù)temp_seriesseries.dropna()# 去除 NaNresultadfuller(temp_series)print(fADF 統(tǒng)計量:{result[0]:.4f})print(fp 值:{result[1]:.4f})print(臨界值:)forkey,valueinresult[4].items():print(f{key}:{value:.4f})ifresult[1]0.05:print(結(jié)論: ★★★ 拒絕原假設(shè)序列是平穩(wěn)的)else:print(結(jié)論: 接受原假設(shè)序列是非平穩(wěn)的 (需要進一步差分))# 對差分后的序列進行檢驗adf_test(df[diff_log_export])# # Step 4: 模型識別與定階 (ACF/PACF)# print( Step 4: 繪制 ACF 和 PACF 圖...)fig,axplt.subplots(2,1,figsize(12,8))# 去除 NaN 值后繪圖plot_acf(df[diff_log_export].dropna(),lags40,axax[0])ax[0].set_title(自相關(guān)函數(shù) (ACF))plot_pacf(df[diff_log_export].dropna(),lags40,axax[1])ax[1].set_title(偏自相關(guān)函數(shù) (PACF))plt.tight_layout()plt.show()# # Step 5: 建立 SARIMA 模型# print( Step 5: 擬合 SARIMA(2,1,1)x(2,1,1,12) 模型...)# 注意這里我們使用原始的 log 數(shù)據(jù)并在 order 中指定 d1# 在 seasonal_order 中指定 D1模型會自動處理差分。# order (p, d, q), seasonal_order (P, D, Q, s)modelsm.tsa.statespace.SARIMAX(df[log_export],order(2,1,1),seasonal_order(2,1,1,12),enforce_stationarityFalse,enforce_invertibilityFalse)resultsmodel.fit(dispFalse)# dispFalse 不打印收斂過程print(模型擬合完成。)print(results.summary().tables[1])# 打印核心參數(shù)表# # Step 6: 模型檢驗 (殘差分析)# print( Step 6: 模型殘差診斷...)# 1. 繪制診斷圖results.plot_diagnostics(figsize(15,12))plt.show()# 2. Ljung-Box 白噪聲檢驗lb_testacorr_ljungbox(results.resid,lags[6,12,24],return_dfTrue)print(Ljung-Box 檢驗結(jié)果 (p值 0.05 則為白噪聲):)print(lb_test)# # Step 7: 模型預(yù)測# print( Step 7: 預(yù)測未來 10 個月...)forecast_steps10pred_ucresults.get_forecast(stepsforecast_steps)# 獲取預(yù)測的均值和置信區(qū)間 (對數(shù)尺度)pred_mean_logpred_uc.predicted_mean pred_ci_logpred_uc.conf_int()# *** 逆變換從對數(shù)還原為原始值 (exp) ***pred_meannp.exp(pred_mean_log)pred_cinp.exp(pred_ci_log)# 生成未來日期索引last_datedf.index[-1]future_datespd.date_range(startlast_datepd.DateOffset(months1),periodsforecast_steps,freqM)pred_mean.indexfuture_dates pred_ci.indexfuture_dates# 打印預(yù)測數(shù)值forecast_dfpd.DataFrame({預(yù)測值:pred_mean,下限 (95%):pred_ci.iloc[:,0],上限 (95%):pred_ci.iloc[:,1]})print(forecast_df)# --- 最終可視化 ---plt.figure(figsize(14,7))# 繪制歷史數(shù)據(jù) (顯示最后 50 個點以便觀察細節(jié))plt.plot(df.index[-50:],df[export][-50:],label歷史觀察值 (最近50個月),colorblack)# 繪制預(yù)測值plt.plot(pred_mean.index,pred_mean,label未來預(yù)測值,colorred,markero,linestyle--)# 繪制置信區(qū)間plt.fill_between(pred_ci.index,pred_ci.iloc[:,0],pred_ci.iloc[:,1],colorpink,alpha0.3,label95% 置信區(qū)間)plt.title(SARIMA 模型未來 10 個月出口額預(yù)測,fontsize16)plt.xlabel(時間)plt.ylabel(出口額)plt.legend(locupper left)plt.grid(True)plt.show()print( 分析結(jié)束。)總結(jié)時間序列分析是一套嚴密的差分方程與隨機過程理論。從平穩(wěn)性的E(Xt)μE(X_t)muE(Xt?)μ到 SARIMA 復(fù)雜的算子多項式Φ(L)ΦS(LS)Phi(L)Phi_S(L^S)Φ(L)ΦS?(LS)每一步都有嚴格的數(shù)學(xué)定義。掌握這些公式能幫助我們更深刻地理解數(shù)據(jù)為何需要差分、參數(shù)代表什么物理意義以及如何科學(xué)地評估模型的有效性。