网站开发 数字证书wordpress 短链接插件

网站开发 数字证书,wordpress 短链接插件,宏大建设集团网站,计算机应用技术培训班时间序列分析#xff1a;从随机过程理论到 SARIMA 模型构建 一、 时间序列的数学定义与性质 1.1 随机序列与观察值 在数学上#xff0c;时间序列不仅仅是一组数字#xff0c;它是一个随机过程#xff08;Stochastic Process#xff09;。 设 TTT 为一个时间索引集合#…时间序列分析从随机过程理论到 SARIMA 模型构建一、 时间序列的数学定义与性质1.1 随机序列与观察值在数学上时间序列不仅仅是一组数字它是一个随机过程Stochastic Process。设TTT为一个时间索引集合通常为整数集Z\mathbb{Z}Z我们称随机变量族{Xt:t∈T}\{X_t : t \in T\}{Xt​:t∈T}为一个时间序列。我们在现实中记录到的具体数据如某年的GDP是这个随机过程的一个样本实现Realization记为{xt:t1,2,…,n}\{x_t : t1, 2, \dots, n\}{xt​:t1,2,…,n}。importpandasaspdimportnumpyasnpimportmatplotlib.pyplotaspltimportseabornassns# 设置绘图风格sns.set(stylewhitegrid)plt.rcParams[font.sans-serif][SimHei]# 用来正常显示中文标签plt.rcParams[axes.unicode_minus]False# 用来正常显示负号# 读取数据# 假设 Excel 文件名为 export_datas.xlsx时间列为 date数据列为 exportdfpd.read_excel(export_datas.xlsx)# 创建时间索引 (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(我国出口额月度数据 (2001-2016),fontsize14)plt.xlabel(年份)plt.ylabel(出口额)plt.legend()plt.show()1.2 平稳性 (Stationarity)平稳性是时间序列建模的基石。如果一个随机过程的统计规律不随时间推移而改变则称其为平稳的。1. 严平稳 (Strict Stationarity)若对任意的t1,t2,…,tk∈Tt_1, t_2, \dots, t_k \in Tt1​,t2​,…,tk​∈T和任意整数τ\tauτ随机变量组(Xt1,…,Xtk)(X_{t_1}, \dots, X_{t_k})(Xt1​​,…,Xtk​​)与(Xt1τ,…,Xtkτ)(X_{t_1\tau}, \dots, X_{t_k\tau})(Xt1​τ​,…,Xtk​τ​)具有相同的联合概率分布即Ft1,…,tk(x1,…,xk)Ft1τ,…,tkτ(x1,…,xk) F_{t_1, \dots, t_k}(x_1, \dots, x_k) F_{t_1\tau, \dots, t_k\tau}(x_1, \dots, x_k)Ft1​,…,tk​​(x1​,…,xk​)Ft1​τ,…,tk​τ​(x1​,…,xk​)则称该序列为严平稳序列。条件极其苛刻实际应用中很少满足。2. 宽平稳 (Weak Stationarity)在实际建模中我们通常使用宽平稳或称二阶平稳。一个时间序列{Xt}\{X_t\}{Xt​}称为宽平稳的若满足以下三个条件均值恒定对于任意ttt期望存在且为常数。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自协方差仅依赖于时间间隔对于任意t,st, st,s协方差仅与时间差ks−tk s-tks−t有关。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​二、 滞后算子与差分为了简化模型的数学表达引入**滞后算子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​常数ccc在滞后算子作用下不变LccL c cLcc。2.2 差分运算差分是处理非平稳序列去除趋势的主要手段。一阶差分Δ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步差分季节差分Δ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阶差分连续进行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 \binom{d}{j} X_{t-j}ΔdXt​j0∑d​(−1)j(jd​)Xt−j​# 1. 对数变换 (稳定方差)df[log_export]np.log(df[export])# 2. 一阶差分 (消除趋势)# 注意差分后会产生 NaN 值需要去除df[diff_log_export]df[log_export].diff().dropna()# 绘制处理后的数据fig,axplt.subplots(2,1,figsize(12,10))# 对数数据图ax[0].plot(df.index,df[log_export],colorgreen)ax[0].set_title(对数变换后的序列)# 差分数据图ax[1].plot(df.index,df[diff_log_export],colororange)ax[1].set_title(对数一阶差分后的序列)plt.tight_layout()plt.show()平稳性检验fromstatsmodels.tsa.stattoolsimportadfullerdefadf_test(series): 封装 ADF 检验函数输出清晰的结果 print(ADF 检验结果:)resultadfuller(series.dropna())print(f统计量 (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(结论: 拒绝原假设序列是平稳的)else:print(结论: 接受原假设序列是非平稳的)# 对差分后的序列进行检验adf_test(df[diff_log_export])三、 随机过程ARIMA 模型是由以下几个基础过程构建而成的。3.1 白噪声过程 (White Noise)白噪声是时间序列建模中的基本“原子”通常记为{εt}\{\varepsilon_t\}{εt​}。若满足以下性质则称为纯随机过程或白噪声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 \neq 0Cov(εt​,εtk​)0,∀k03.2 自回归过程 AR§ppp阶自回归模型 (AutoRegressive model of orderppp)描述了当前值与过去ppp个历史值之间的线性关系。数学表达式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为常数项ϕi\phi_iϕi​为自回归系数εt\varepsilon_tεt​为白噪声。引入滞后算子多项式Φ(L)\Phi(L)Φ(L)令Φ(L)1−ϕ1L−ϕ2L2−⋯−ϕpLp\Phi(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​平稳性条件AR§ 序列平稳的充要条件是特征方程Φ(z)0\Phi(z) 0Φ(z)0的所有根zzz都在单位圆外即∣z∣1,for all z such that 1−∑i1pϕizi0 |z| 1, \quad \text{for all } z \text{ 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)描述了当前值与过去qqq个随机扰动项误差之间的线性关系。数学表达式Xtμεt−θ1εt−1−θ2εt−2−⋯−θqεt−q X_t \mu \varepsilon_t - \theta_1 \varepsilon_{t-1} - \theta_2 \varepsilon_{t-2} - \dots - \theta_q \varepsilon_{t-q}Xt​μεt​−θ1​εt−1​−θ2​εt−2​−⋯−θq​εt−q​(注根据 Box-Jenkins 习惯MA 参数前常取负号也有教材取正号本质只需调整参数符号即可)引入滞后算子多项式Θ(L)\Theta(L)Θ(L)令Θ(L)1−θ1L−θ2L2−⋯−θqLq\Theta(L) 1 - \theta_1 L - \theta_2 L^2 - \dots - \theta_q L^qΘ(L)1−θ1​L−θ2​L2−⋯−θq​Lq则模型可简写为XtμΘ(L)εt X_t \mu \Theta(L)\varepsilon_tXt​μΘ(L)εt​性质MA(q) 模型总是平稳的。但为了保证模型的可逆性 (Invertibility)即可以用过去的观察值表示当前的误差要求方程Θ(z)0\Theta(z) 0Θ(z)0的所有根都在单位圆外。3.4 自回归移动平均模型 ARMA(p, q)将 AR 和 MA 结合得到 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 \theta_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))# 自相关图 (ACF)plot_acf(df[diff_log_export].dropna(),lags40,axax[0])ax[0].set_title(自相关函数 (ACF))# 偏自相关图 (PACF)plot_pacf(df[diff_log_export].dropna(),lags40,axax[1])ax[1].set_title(偏自相关函数 (PACF))plt.tight_layout()plt.show()# 提示也可以使用 pmdarima 库进行自动定阶# import pmdarima as pm# model pm.auto_arima(df[log_export], seasonalTrue, m12)四、 ARIMA 与 SARIMA非平稳与季节性建模4.1 ARIMA(p, d, q) 模型现实数据往往是非平稳的含有趋势。如果一个时间序列{Yt}\{Y_t\}{Yt​}经过ddd次差分后变为平稳的 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)εt\Phi(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\right) (1 - L)^d Y_t c \left(1 - \sum_{j1}^q \theta_j L^j\right) \varepsilon_t(1−i1∑p​ϕi​Li)(1−L)dYt​c(1−j1∑q​θj​Lj)εt​4.2 SARIMA 模型季节性的引入对于具有周期性如月度周期S12S12S12季度周期S4S4S4的数据普通的 ARIMA 无法捕捉季节效应。我们需要引入季节性自回归移动平均模型 (SARIMA)。SARIMA 模型的结构由非季节部分(p,d,q)(p, d, q)(p,d,q)和季节部分(P,D,Q)S(P, D, Q)_S(P,D,Q)S​组成记为ARIMA(p,d,q)×(P,D,Q)S ARIMA(p, d, q) \times (P, D, Q)_SARIMA(p,d,q)×(P,D,Q)S​1. 季节算子定义季节性 AR 多项式(阶数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季节性 MA 多项式(阶数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季节性差分算子(阶数DDD)ΔSD(1−LS)D \Delta_S^D (1 - L^S)^DΔSD​(1−LS)D2. SARIMA 通用数学表达式Φ(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. 公式的详细结构解析为了完整展示该模型我们将上述各项代入并展开(1−∑i1pϕiLi)⏟非季节 AR(1−∑j1PΦjLjS)⏟季节 AR(1−L)d⏟差分(1−LS)D⏟季节差分Yt(1−∑k1qθkLk)⏟非季节 MA(1−∑m1QΘmLmS)⏟季节 MAεt \underbrace{\left(1 - \sum_{i1}^p \phi_i L^i\right)}_{\text{非季节 AR}} \underbrace{\left(1 - \sum_{j1}^P \Phi_j L^{jS}\right)}_{\text{季节 AR}} \underbrace{(1-L)^d}_{\text{差分}} \underbrace{(1-L^S)^D}_{\text{季节差分}} Y_t \underbrace{\left(1 - \sum_{k1}^q \theta_k L^k\right)}_{\text{非季节 MA}} \underbrace{\left(1 - \sum_{m1}^Q \Theta_m L^{mS}\right)}_{\text{季节 MA}} \varepsilon_t非季节AR(1−i1∑p​ϕi​Li)​​季节AR(1−j1∑P​Φj​LjS)​​差分(1−L)d​​季节差分(1−LS)D​​Yt​非季节MA(1−k1∑q​θk​Lk)​​季节MA(1−m1∑Q​Θm​LmS)​​εt​这是时间序列分析中最为复杂也最为通用的线性模型形式。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()# 打印模型详细统计摘要print(results.summary())五、 参数估计与模型选择5.1 极大似然估计 (MLE)在确定了模型阶数(p,d,q)(p, d, q)(p,d,q)和(P,D,Q)(P, D, Q)(P,D,Q)后需要估计参数向量β(ϕ,θ,Φ,Θ,σ2)\beta (\phi, \theta, \Phi, \Theta, \sigma^2)β(ϕ,θ,Φ,Θ,σ2)。通常假设残差服从正态分布εt∼N(0,σ2)\varepsilon_t \sim N(0, \sigma^2)εt​∼N(0,σ2)使用对数似然函数ln⁡L(β∣Y)≈−n2ln⁡(2π)−n2ln⁡(σ2)−12σ2∑t1nεt2(β) \ln L(\beta | Y) \approx -\frac{n}{2}\ln(2\pi) - \frac{n}{2}\ln(\sigma^2) - \frac{1}{2\sigma^2} \sum_{t1}^n \varepsilon_t^2(\beta)lnL(β∣Y)≈−2n​ln(2π)−2n​ln(σ2)−2σ21​t1∑n​εt2​(β)通过数值优化算法如 BFGS 或 Newton-Raphson最大化该函数以求解参数。5.2 信息准则 (Information Criteria)为了在模型拟合优度与模型复杂度之间取得平衡避免过拟合我们使用信息准则进行模型选择。1. 赤池信息准则 (AIC)AIC−2ln⁡(L^)2k AIC -2 \ln(\hat{L}) 2kAIC−2ln(L^)2k其中L^\hat{L}L^是似然函数的最大值kkk是模型中被估计参数的总数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 对参数个数的惩罚力度比 AIC 更大当n≥8n \geq 8n≥8时倾向于选择更简单的模型。判别规则AIC 或 BIC 值越小模型越优。六、 残差检验建模的最后一步是验证模型假设即检验残差序列{ε^t}\{\hat{\varepsilon}_t\}{ε^t​}是否为白噪声。如果残差中仍包含自相关性说明模型未能提取全部信息。Ljung-Box Q 检验原假设H0H_0H0​残差序列不存在自相关即ρ1ρ2⋯ρm0\rho_1 \rho_2 \dots \rho_m 0ρ1​ρ2​⋯ρm​0。检验统计量QQQ定义为Qn(n2)∑k1mρ^k2n−k Q n(n2) \sum_{k1}^m \frac{\hat{\rho}_k^2}{n-k}Qn(n2)k1∑m​n−kρ^​k2​​其中nnn是样本量。ρ^k\hat{\rho}_kρ^​k​是残差的kkk阶样本自相关系数ρ^k∑tk1nε^tε^t−k∑t1nε^t2 \hat{\rho}_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是检验的最大滞后阶数。统计量QQQ近似服从自由度为m−(pqPQ)m - (pqPQ)m−(pqPQ)的卡方分布Q∼χ2(m−params) Q \sim \chi^2(m - \text{params})Q∼χ2(m−params)若计算出的ppp值小于显著性水平α\alphaα如 0.05则拒绝原假设认为模型不合格反之则认为模型有效。fromstatsmodels.stats.diagnosticimportacorr_ljungbox# 1. 绘制残差诊断图 (包括直方图、QQ图、相关图)results.plot_diagnostics(figsize(15,12))plt.show()# 2. Ljung-Box 白噪声检验# 检验残差是否存在自相关选取滞后 6 期和 12 期lb_testacorr_ljungbox(results.resid,lags[6,12],return_dfTrue)print(Ljung-Box 检验结果:)print(lb_test)# 如果 lb_pvalue 0.05则说明残差为白噪声模型通过检验模型预测# 预测未来 10 个月forecast_steps10pred_ucresults.get_forecast(stepsforecast_steps)# 获取预测的均值和置信区间pred_mean_logpred_uc.predicted_mean pred_ci_logpred_uc.conf_int()# *** 关键步骤逆变换 (从对数还原为原始值) ***pred_meannp.exp(pred_mean_log)pred_cinp.exp(pred_ci_log)# 创建未来时间索引future_datespd.date_range(startdf.index[-1]pd.DateOffset(months1),periodsforecast_steps,freqM)pred_mean.indexfuture_dates pred_ci.indexfuture_dates# 打印预测结果forecast_dfpd.DataFrame({预测值:pred_mean,下限 (95%):pred_ci.iloc[:,0],上限 (95%):pred_ci.iloc[:,1]})print(forecast_df)# 绘图展示plt.figure(figsize(14,7))# 绘制历史数据plt.plot(df.index,df[export],label历史观察值,colorblack)# 绘制预测值plt.plot(pred_mean.index,pred_mean,label未来预测值,colorred,markero)# 绘制置信区间plt.fill_between(pred_ci.index,pred_ci.iloc[:,0],pred_ci.iloc[:,1],colorpink,alpha0.3,label95% 置信区间)plt.title(SARIMA 模型未来 10 个月出口额预测,fontsize16)plt.legend(locupper left)plt.grid(True)plt.show()七、 案例数学模型假设我们最终确定的模型为ARIMA(2,1,1)×(2,1,1)12ARIMA(2,1,1) \times (2,1,1)_{12}ARIMA(2,1,1)×(2,1,1)12​且对原始数据YtY_tYt​进行了对数变换ytln⁡(Yt)y_t \ln(Y_t)yt​ln(Yt​)。这意味着我们需要对yty_tyt​进行一次普通差分(1−L)(1-L)(1−L)和一次周期为 12 的季节差分(1−L12)(1-L^{12})(1−L12)。该模型的具体数学方程展开如下(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 - \theta_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​)对数变换解决异方差。(1−L)(1−L12)(1-L)(1-L^{12})(1−L)(1−L12)双重差分消除长期趋势和季节性趋势。(1−ϕ1L−ϕ2L2)(1 - \phi_1 L - \phi_2 L^2)(1−ϕ1​L−ϕ2​L2)当前的差分值与过去 1 个月和 2 个月的差分值相关。(1−Φ1L12−Φ2L24)(1 - \Phi_1 L^{12} - \Phi_2 L^{24})(1−Φ1​L12−Φ2​L24)当前的季节性变化与去年的同月、前年的同月的变化相关。右边(1−θ1L)(1 - \theta_1 L)(1−θ1​L)包含 1 阶非季节性移动平均成分短期震荡。(1−Θ1L12)(1 - \Theta_1 L^{12})(1−Θ1​L12)包含 1 阶季节性移动平均成分季节性震荡。εt\varepsilon_tεt​剩余的纯随机波动白噪声。通过求解上述方程的期望即可得到对未来时刻YnhY_{nh}Ynh​的最优预测值。八、完整代码# -*- coding: utf-8 -*- 时间序列分析完整实战代码 (SARIMA模型) 对应博客章节七、 Python 实战我国出口额月度数据预测 importpandasaspdimportnumpyasnpimportmatplotlib.pyplotaspltimportseabornassnsimportstatsmodels.apiassmfromstatsmodels.tsa.stattoolsimportadfullerfromstatsmodels.graphics.tsaplotsimportplot_acf,plot_pacffromstatsmodels.stats.diagnosticimportacorr_ljungboximportwarnings# 忽略一些不必要的警告warnings.filterwarnings(ignore)# # 0. 全局设置# # 设置绘图风格sns.set(stylewhitegrid)# 设置字体以支持中文显示# Windows用户通常使用 SimHeiMac用户使用 Arial Unicode MSplt.rcParams[font.sans-serif][SimHei]plt.rcParams[axes.unicode_minus]False# 开关是否使用模拟数据# 如果你有真实的 export_datas.xlsx 文件请将此处改为 FalseUSE_MOCK_DATATrue# # Step 1: 读取数据与初步观察# print( Step 1: 数据准备中...)ifUSE_MOCK_DATA:# --- 生成模拟数据 (为了演示代码可运行) ---datespd.date_range(start2001-01-01,periods188,freqM)# 模拟趋势 季节性 随机噪声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# 保证为正数dfpd.DataFrame({date:dates,export:data_values})df.set_index(date,inplaceTrue)print(提示正在使用模拟数据运行...)else:# --- 读取真实数据 ---try:dfpd.read_excel(export_datas.xlsx)# 假设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 设置为 True 以测试代码。)exit()# 绘制原始时序图plt.figure(figsize(12,6))plt.plot(df.index,df[export],label出口额(亿美元),colornavy)plt.title(我国出口额月度数据 (2001-2016),fontsize14)plt.xlabel(年份)plt.ylabel(出口额)plt.legend()plt.show()# # Step 2: 数据平稳化处理# print(\n Step 2: 进行对数变换与差分处理...)# 1. 对数变换 (稳定方差)df[log_export]np.log(df[export])# 2. 一阶差分 (消除趋势)df[diff_log_export]df[log_export].diff()# 绘制处理后的数据对比图fig,axplt.subplots(2,1,figsize(12,10))ax[0].plot(df.index,df[log_export],colorgreen)ax[0].set_title(对数变换后的序列)ax[0].grid(True)ax[1].plot(df.index,df[diff_log_export],colororange)ax[1].set_title(对数一阶差分后的序列)ax[1].grid(True)plt.tight_layout()plt.show()# # Step 3: 平稳性检验 (ADF Test)# print(\n Step 3: 进行 ADF 平稳性检验...)defadf_test(series):ADF 检验封装函数temp_seriesseries.dropna()# 去除 NaNresultadfuller(temp_series)print(fADF 统计量:{result[0]:.4f})print(fp 值:{result[1]:.4f})print(临界值:)forkey,valueinresult[4].items():print(f{key}:{value:.4f})ifresult[1]0.05:print(结论: ★★★ 拒绝原假设序列是平稳的)else:print(结论: 接受原假设序列是非平稳的 (需要进一步差分))# 对差分后的序列进行检验adf_test(df[diff_log_export])# # Step 4: 模型识别与定阶 (ACF/PACF)# print(\n 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(自相关函数 (ACF))plot_pacf(df[diff_log_export].dropna(),lags40,axax[1])ax[1].set_title(偏自相关函数 (PACF))plt.tight_layout()plt.show()# # Step 5: 建立 SARIMA 模型# print(\n Step 5: 拟合 SARIMA(2,1,1)x(2,1,1,12) 模型...)# 注意这里我们使用原始的 log 数据并在 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])# 打印核心参数表# # Step 6: 模型检验 (残差分析)# print(\n 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 检验结果 (p值 0.05 则为白噪声):)print(lb_test)# # Step 7: 模型预测# print(\n Step 7: 预测未来 10 个月...)forecast_steps10pred_ucresults.get_forecast(stepsforecast_steps)# 获取预测的均值和置信区间 (对数尺度)pred_mean_logpred_uc.predicted_mean pred_ci_logpred_uc.conf_int()# *** 逆变换从对数还原为原始值 (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# 打印预测数值forecast_dfpd.DataFrame({预测值:pred_mean,下限 (95%):pred_ci.iloc[:,0],上限 (95%):pred_ci.iloc[:,1]})print(forecast_df)# --- 最终可视化 ---plt.figure(figsize(14,7))# 绘制历史数据 (显示最后 50 个点以便观察细节)plt.plot(df.index[-50:],df[export][-50:],label历史观察值 (最近50个月),colorblack)# 绘制预测值plt.plot(pred_mean.index,pred_mean,label未来预测值,colorred,markero,linestyle--)# 绘制置信区间plt.fill_between(pred_ci.index,pred_ci.iloc[:,0],pred_ci.iloc[:,1],colorpink,alpha0.3,label95% 置信区间)plt.title(SARIMA 模型未来 10 个月出口额预测,fontsize16)plt.xlabel(时间)plt.ylabel(出口额)plt.legend(locupper left)plt.grid(True)plt.show()print(\n 分析结束。)总结时间序列分析是一套严密的差分方程与随机过程理论。从平稳性的E(Xt)μE(X_t)\muE(Xt​)μ到 SARIMA 复杂的算子多项式Φ(L)ΦS(LS)\Phi(L)\Phi_S(L^S)Φ(L)ΦS​(LS)每一步都有严格的数学定义。掌握这些公式能帮助我们更深刻地理解数据为何需要差分、参数代表什么物理意义以及如何科学地评估模型的有效性。