建议：数据的预处理（第一部分）使用Python，模型的定阶及预测（第二部分）使用SAS。
一、时间序列的预处理
1.时间序列的平稳性检验及平稳化处理(Python)
首先，对一个时间序列进行ARIMA（p,d,q）模型建模，需要判断样本的平稳性，如果不是平稳序列，需要利用一次或多次差分将其转化为平稳序列。
对于平稳性检验，python做的较好，示例如下（工具ipython jupyter notebook）：
1)使用Pandas进行时间序列处理
#所需头文件
import pandas as pd
import numpy as np
import matplotlib.pylab as plt
%matplotlib inline
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6
2）从文件读取数据
# dateparse('1962-01')
dateparse = lambda dates: pd.datetime.strptime(dates, '%Y-%m')
#这样得到的才是时间序列(得到的才是一个 TS object)
data = \
pd.read_csv('AirPassengers.csv', parse_dates='Month', index_col='Month',date_parser=dateparse)
print data.head()
3）Optional：Convert the column into a Series object to prevent referring to columns names every time when using the TS.
#convert to time series:
ts = data['#Passengers']
ts.head(10)
4)使用移动平均和移动标准差以及Dickey Fuller Test检查序列的平稳性
平稳性的三个标准：
- constant mean
- constant variance
- an autocovariance that does not depend on time.
I)平稳性感性识别
plt.plot(ts,color='gray',linestyle='--',marker='o')
plt.grid()
II)定义平稳性测试函数
#平稳性测试函数
from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries):
#Determing rolling statistics
#计算移动平均，窗口大小为12
rolmean = pd.rolling_mean(timeseries, window=12)
#计算移动标准差，窗口大小为12
rolstd = pd.rolling_std(timeseries, window=12)
#Plot rolling statistics:
orig = plt.plot(timeseries, color='blue',label='Original')
mean = plt.plot(rolmean, color='red', label='Rolling Mean')
std = plt.plot(rolstd, color='black', label = 'Rolling Std')
plt.legend(loc='best')
plt.title('Rolling Mean & Standard Deviation')
plt.show(block=False)
#Perform Dickey-Fuller test:
print 'Results of Dickey-Fuller Test:'
dftest = adfuller(timeseries, autolag='AIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
dfoutput['Critical Value (%s)'%key] = value
print dfoutput
III)测试时间序列的平稳性
test_stationarity(ts)
输出结果
Results of Dickey-Fuller Test:
Test Statistic 0.815369
p-value 0.991880
#Lags Used 13.000000
Number of Observations Used 130.000000
Critical Value (5%) -2.884042
Critical Value (1%) -3.481682
Critical Value (10%) -2.578770
dtype: float64
1）我们需要将Test Statistic值（即Pdf值）与三个Critical Value值进行比较（有符号数比较），小于某一个值即可说明有1-α的把握认为序列平稳，α=5%/1%/10% 。
2）p-value越小越可以拒绝原假设，原假设是序列不平稳。
由输出结果可知，ts不平稳。
5）平稳化处理
有许多种方法可以使得序列平稳化，这里主要介绍差分方法。
因为序列有明显一直增长的趋势，所以这里先使用log减少这种趋势：
ts_log = np.log(ts)
ts_log_diff = ts_log.diff()# 可使用ts_log.diff(2)进行二阶差分
plt.plot(ts_log_diff)
差分后检验平稳性：
ts_log_diff.dropna(inplace=True)
test_stationarity(ts_log_diff)
此时，有输出结果可知，拥有90%的把握认为序列是平稳的。
什么是差分？一阶差分指原序列值相距一期的两个序列值之间的减法运算；k阶差分就是相距k期的两个序列值之间相减。如果一个时间序列经过差分运算后具有平稳性，则该序列为差分平稳序列，可以使用ARIMA模型进行分析。
2.零均值检验(Python)
一个零均值平稳时间序列的自相关和偏自相关函数，要么拖尾，要么截尾。如果零值化的时序既不拖尾，也不截尾，而是呈现出缓慢衰减或者周期性衰减，则认为可能存在趋势或周期性，应视为非平稳。
对样本的均值是否为零进行检验。如果显著不为零，就需要对序列进行零均值化。
零均值的检验假设：H0:μ=EXt=0。
设时间序列的均值为μ，标准差为ρ。如果满足-2ρ<=μ<=2ρ,那么接受原假设，即认为均值显著为0.
疑惑：网上大多例子中并没有进行零均值检验，感觉即使均值不为零也可建模。
I)定义零均值测试函数
#定义零均值测试函数
def test_mean(timeseries):
mu=timeseries.mean()
std=timeseries.std()
print u'均值 =',mu
print u'方差 =',std
plt.axhline(y=-2*std,linestyle='--',color='gray')
plt.axhline(y=mu,linestyle='-',color='red')
plt.axhline(y=2*std,linestyle='--',color='gray')
plt.title('Autocorrelation Function')
if mu >= -2*std and mu <= 2*std:
print '零均值检验结果：均值显著为0.'
else:
print '零均值检验结果：均值不显著为0.'
II)进行零均值检验
#进行零均值检验
test_mean(ts_log_diff)
3.白噪声检验(Python)
对于纯随机序列，又称白噪声序列，序列的各项数值之间没有任何相关关系，序列在进行完全无序的随机波动，可以终止对该序列的分析。白噪声序列是没有信息可提取的平稳序列。
白噪声检验原假设：时间序列是白噪声。
使用python示例如下：
#白噪声检验
from statsmodels.stats.diagnostic import acorr_ljungbox
#返回统计量和p值
noiseRes = acorr_ljungbox(ts_log_diff, lags=1)
print u'一阶差分序列的白噪声检验结果为：'
print 'stat | p-value'
for x in noiseRes:
print x,'|',
输出结果如下：
一阶差分序列的白噪声检验结果为：
stat | p-value
[ 5.82633006] | [ 0.01578803] |
P值小于0.05，拒绝原假设，所以一阶差分后的序列为平稳非白噪声序列(95%的把握)。
对于平稳非白噪声序列，它的均值和方差是常数。通常是建立一个线性模型来拟合该序列的发展，借此提取该序列的有用信息。ARIMA（p,d,q）模型是最常用的平稳序列拟合模型(加入I主要是将非平稳序列转换为平稳序列，d是差分阶数)。
二、平稳时间序列建模
4.模型初步识别(SAS & python)
根据序列的自相关函数和偏自相关函数的特征可以初步判断模型类型，如下表：
| 自相关函数（ACF） | 偏自相关函数（PACF） | 选择模型 |
|---|---|---|
| 拖尾 | p阶截尾 | AR(p) |
| q阶截尾 | 拖尾 | MA(q) |
| p阶拖尾 | q阶拖尾 | ARMA(p，q) |
ARIMA 是 ARMA 算法的扩展版，用法类似 。
但是，如果样本的自相关函数ρ_k和偏自相关函数φ_{kk}没有以上特点，随着k的增加衰减缓慢或呈现周期性衰减，说明该序列仍然是不平稳的，需要进行序列的平稳化处理。
参考资料1：
自回归模型（AR模型）、移动平均模型（MA模型）和自回归移动平均模型（ARMA模型）阶数识别，确定模型阶数p和q值：
- AR模型：某个观测值Xt与其滞后p期的观测值的线性组合再加上随机误差项。
即：Xt = φ1Xt-1 + φ2Xt-2 + …… + φpXt-p + at；- MA模型：某个观测值Xt与先前t-1,t-2,t-q个时刻进入系统的q个随机误差项即at，at-1，……，at-q的线性组合。
即：Xt = at - θ1at-1 - θ2at-2 - …… - θqat-q；- ARMA模型：即观测值不仅与其以前p个时刻的自身观测值有关，而且还与其以前时刻进入系统的q个随机误差存在一定的依存关系。
即Xt= φ1Xt-1 + φ2Xt-2 + …… + φpXt-p + at - θ1at-1 - θ2at-2- …… - θqat-q。
参考资料2：
These can be used to determine the ‘p’ and ‘q’ values as:
- p – The lag value where the PACF chart crosses the upper confidence interval for the first time. If you notice closely, in this case p=2.
- q – The lag value where the ACF chart crosses the upper confidence interval for the first time. If you notice closely, in this case q=2.
参考资料3：
ACF plot is a bar chart of the coefficients of correlation between a time series and lags of itself.
PACF plot is a plot of the partial correlation coefficients between the series and lags of itself.
To find p and q you need to look at ACF and PACF plots. The interpretation of ACF and PACF plots to find p and q are as follows:
- AR (p) model: If ACF plot tails off but PACF plot cut off after p lags
- MA(q) model: If PACF plot tails off but ACF plot cut off after q lags
- ARMA(p,q) model: If both ACF and PACF plot tail off, you can choose different combinations of p and q , smaller p and q are tried.
- ARIMA(p,d,q) model: If it’s ARMA with d times differencing to make time series stationary.
Use AIC and BIC to find the most appropriate model. Lower values of AIC and BIC are desirable.
Tails off means slow decaying of the plot, i.e. plot has significant spikes at higher lags too. Cut off means the bar is significant at lag p and not significant at any higher order lags.
当然，以上只是对p和q的初步估计，确定p和q在第5步。
ACF和PACF的刻画示例如下：
1）在此，先介绍使用SAS进行模型初步识别：
首先导入数据，根据第一部分的数据预处理方法对数据进行相同的预处理：
proc import datafile= '/folders/myfolders/data/AirPassengers.csv' dbms = csv out = simple_air;
run;
data air;
set simple_air;
rename _Passengers=P;
plog = log(_Passengers);/*进行log处理*/
run;
/* 画出序列趋势图 */
proc sgplot data=air;
title "Air Passengers";
series x=Month y=plog / legendlabel = 'Passengers'
markers lineattrs = (color = red pattern = solid thickness = 0.5)
markerattrs=(color = red symbol=circlefilled);
run;
/* 模型初步识别 */
proc arima data=air;
identify var=plog(1) nlag=30;/*plog(1)即是进行一阶差分，nlag表示滞后的阶数*/
run;
2）使用python刻画序列的自相关函数和偏自相关函数如下:
#自相关图
from statsmodels.graphics.tsaplots import plot_acf
plot_acf(ts_log_diff).show()
plt.show()
from statsmodels.graphics.tsaplots import plot_pacf
#偏自相关图
plot_pacf(ts_log_diff).show()
plt.show()
或者
#ACF and PACF plots:
from statsmodels.tsa.stattools import acf, pacf
lag_acf = acf(ts_log_diff, nlags=20)
lag_pacf = pacf(ts_log_diff, nlags=20, method='ols')
#Plot ACF:
plt.subplot(121)
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Autocorrelation Function')
#Plot PACF:
plt.subplot(122)
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Partial Autocorrelation Function')
plt.tight_layout()
5.模型定阶（SAS & python）
确定模型的类型后，下一步就是确定模型的阶数。
对于AR(q)和MA(q)这两个简单模型，利用自相关函数和偏自相关函数，在确定模型类型的同时，就可以初步判断阶数，但仍然无法确定其精确的阶数。
对于复杂的ARMA(p,q)模型，难以用自相关函数和偏自相关函数确定阶数，可以使用Akaike(1976)提出的BIC准则，也可以使用AIC准则（都是值越小越好）。
此处介绍使用SAS和python两种方法，推荐使用SAS,其处理速度更快。
1）使用SAS定阶：
/* 模型定阶/
proc arima data=air;
identify var=plog(1) nlag=30 minic p=(0:14) q=(0:14);
/* 还可以添加选项minic, esacf, scan */
run;
在输出结果的“Minimum Information Criterion”中，可以看到“Minimum Table Value：BIC(p,q)”,其中BIC(p,q)中的值p,q即为所需的AR和MA阶数。
2）使用python定阶(效率较低，不建议使用):
# statsmodels.tsa.stattools.arma_order_select_ic(y, max_ar=4, max_ma=2, ic='bic', trend='c', model_kw={}, fit_kw={})
from statsmodels.tsa.stattools import arma_order_select_ic
res = arma_order_select_ic(ts_log_diff,3,3,ic=['aic','bic'],trend='nc')
# res.aic_min_order
res.bic_min_order
输出(p,q)值即为所求的AR和MA阶数。
6.模型参数估计（SAS）
在模型按照确定的阶数后，接下来就是确定模型的参数了。参数估计主要有三种方法：矩估计、最小二乘法和极大似然法。这三种方法各有利弊，为达到最佳的模型拟合效果，可能需要对这三种方法进行尝试，最后根据样本数据选取最优估计方法。
使用SAS进行模型参数估计示例：
/*使用极大似然估计法进行参数估计*/
proc arima data=air;
identify var=plog(1) nlag=30;
estimate p=13 q=0 method=ml;
/* 还可以添加选项method=ml（极大似然）、uls（非条件最小二乘法）、cls（最小二乘法） */
run;
7.参数显著性检验(SAS)
在建立模型之后，需要进行模型参数的显著性检验（即检验模型中每个未知参数是否显著为零）和残差的白噪声检验（原假设为残差是白噪）。经过参数的显著性检验后，将参数不显著的项从模型中删除，并重新建模，直至得到一个更简单的模型。
原则如下：
(1)在第五步SAS输出结果的“Maximum Likelihood Estimation（或Conditional Least Squares Estimation或Unconditional Least Squares Estimation）”中，查看所有被估计参数的检验P值（“Approx Pr > |t|”项）是否小于0.0001(其实只要p值都小于0.05就可说明显著不为0,当然越小越显著不为0)，小于则拒绝参数为0的原假设，即认为未知参数显著。
(2)在“Autocorrelation Check for Residuals”中，查看P值（Pr > ChiSq）是否都明显大于0.05(显示为“.”则不大于0.05)，都大于0.05则认为残差序列为白噪声序列（接受原假设），并认为模型拟合良好。
使用原则（1）来简化模型；使用原则（2）来判定模型是否拟合良好。两者相辅相成，但是首要保证原则（2）。
在首要满足原则（2）并使用原则（1）后，经过反复调整p,q值以及method（ml/cls/uls），得到如下的模型拟合方式：
/* 模型拟合 */
proc arima data=air;
identify var=plog(1) nlag=30;
estimate p=(8,12) q=(1,3) method=cls;/*p利用Lag8和Lag12,q利用Lag1和Lag3*/
/* 还可以添加选项method=ML（极大似然）、ULS（非条件最小二乘法）、CLS（最小二乘法） */
run;
此时满足原则(2)，但是并非完全满足原则(1)。
8.模型预测（SAS）
目前对平稳时间序列的常用预测方法很多，较常用的是条件期望预测法。
在SAS中进行模型预测如下：
/* 模型预测 */
proc arima data=air plots(only)=forecast(forecast) /*plots=all*/;
identify var=plog(1) nlag=30;
estimate p=(8,12) q=(1,3) method=cls;
/*不要使用id属性，因为id增加会按照最小时间间隔增加（这里是秒，不是按月），
这样画出的预测图上的多步预测值几乎都在挤在一起。*/
forecast lead=12 out=out;/*lead=12表示向前预测12步，并输出到数据out*/
run;
将预测数据和实际数据做对比：
/* 模型预测(获得带id的数据) */
proc arima data=air;
identify var=plog(1) nlag=30;
estimate p=(8,12) q=(1,3) method=cls;
forecast lead=12 id=Month out=outAir;
run;
/* 将实际数据和预测数据放在一张图上做对比 */
proc sgplot data=outAir;
series x=Month y=plog / legendlabel = 'actual'
markers lineattrs = (color=black pattern=dash thickness = 0.5)
markerattrs=(color=black symbol=diamond);
series x=Month y=forecast / legendlabel = 'forecast'
markers lineattrs = (color=red pattern=solid thickness = 0.5)
markerattrs=(color=red symbol=dot);
/* 置信区间为95%的下界 */
/* series x=Month y=l95 / legendlabel = 'l95' */
/* markers lineattrs = (color=green pattern=1 thickness = 0.5) */
/* markerattrs=(color=green symbol=star); */
/* 置信区间为95%的上界 */
/* series x=Month y=u95 / legendlabel = 'u95' */
/* markers lineattrs = (color=blue thickness = 0.5) */
/* markerattrs=(color=blue symbol=circlefield); */
title'Air Passengers';
run;
三、写在后面
9.参考资料
- A Complete Tutorial on Time Series Modeling in R
- Complete guide to create a Time Series Forecast(with Codes in Python)
- 时间序列(ARIMA)模型
- 用SAS操作ARIMA时间序列模型的例子
- SAS时间序列模型预测未来航班数量
- 时间序列ARIMA模型详解：python实现店铺一周销售量预测
10.附录
（1）python数据预处理源程序（使用Jupyter notebook）
#1.序列的平稳性检验及平稳化处理
#所需头文件
import pandas as pd
import numpy as np
import matplotlib.pylab as plt
%matplotlib inline
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6
#导入数据
# dateparse('1962-01')
dateparse = lambda dates: pd.datetime.strptime(dates, '%Y-%m')
#这样得到的才是时间序列(得到的才是一个 TS object)
data = \
pd.read_csv('AirPassengers.csv', parse_dates='Month', index_col='Month',date_parser=dateparse)#这样得到的才是时间序列(得到的才是一个 TS object)
print data.head()
#转化为时间序列对象
ts = data['#Passengers']
ts.head(10)
#画出数据散点图
plt.plot(ts,color='gray',linestyle='--',marker='o')
plt.grid()
#定义平稳性测试函数
from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries):
#Determing rolling statistics
rolmean = pd.rolling_mean(timeseries, window=12)
rolstd = pd.rolling_std(timeseries, window=12)
#Plot rolling statistics:
orig = plt.plot(timeseries, color='blue',label='Original')
mean = plt.plot(rolmean, color='red', label='Rolling Mean')
std = plt.plot(rolstd, color='black', label = 'Rolling Std')
plt.legend(loc='best')
plt.title('Rolling Mean & Standard Deviation')
plt.show(block=False)
#Perform Dickey-Fuller test:
print 'Results of Dickey-Fuller Test:'
dftest = adfuller(timeseries, autolag='AIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
dfoutput['Critical Value (%s)'%key] = value
print dfoutput
#进行平稳性测试
test_stationarity(ts)
#平稳化处理（差分）
ts_log = np.log(ts)
ts_log_diff = ts_log.diff()# 可使用ts_log.diff(2)进行二阶差分
plt.plot(ts_log_diff)
#再次进行平稳性测试
ts_log_diff.dropna(inplace=True)
test_stationarity(ts_log_diff)
# 2.零均值检验
#定义零均值测试函数
def test_mean(timeseries):
mu=timeseries.mean()
std=timeseries.std()
print u'均值 =',mu
print u'方差 =',std
plt.axhline(y=-2*std,linestyle='--',color='gray')
plt.axhline(y=mu,linestyle='-',color='red')
plt.axhline(y=2*std,linestyle='--',color='gray')
plt.title('Autocorrelation Function')
if mu >= -2*std and mu <= 2*std:
print '零均值检验结果：均值显著为0.'
else:
print '零均值检验结果：均值不显著为0.'
#进行零均值检验
test_mean(ts_log_diff)
# 3.白噪声检验
#进行白噪声检验
from statsmodels.stats.diagnostic import acorr_ljungbox
#返回统计量和p值
noiseRes = acorr_ljungbox(ts_log_diff, lags=1)
print u'二阶差分序列的白噪声检验结果为：'
print 'stat | p-value'
for x in noiseRes:
print x,'|',
# 4.模型初步识别
#自相关图
from statsmodels.graphics.tsaplots import plot_acf
plot_acf(ts_log_diff).show()
plt.show()
from statsmodels.graphics.tsaplots import plot_pacf
#偏自相关图
plot_pacf(ts_log_diff).show()
plt.show()
#ACF and PACF plots:
from statsmodels.tsa.stattools import acf, pacf
lag_acf = acf(ts_log_diff, nlags=20)
lag_pacf = pacf(ts_log_diff, nlags=20, method='ols')
#Plot ACF:
plt.subplot(121)
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Autocorrelation Function')
#Plot PACF:
plt.subplot(122)
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Partial Autocorrelation Function')
plt.tight_layout()
# 5.模型定阶
# statsmodels.tsa.stattools.arma_order_select_ic(y, max_ar=4, max_ma=2, ic='bic', trend='c', model_kw={}, fit_kw={})
from statsmodels.tsa.stattools import arma_order_select_ic
res = arma_order_select_ic(ts_log_diff,14,3,ic=['aic','bic'],trend='nc')
# res.aic_min_order
res.bic_min_order
（2）完整SAS源程序
proc import datafile= '/folders/myfolders/data/AirPassengers.csv' dbms = csv out = simple_air;
run;
/* 对数据进行对数处理 */
data air;
set simple_air;
rename _Passengers=P;
plog = log(_Passengers);/*进行log处理*/
run;
/* 画出序列趋势图 */
proc sgplot data=air;
title "Air Passengers";
series x=Month y=plog / legendlabel = 'Passengers'
markers lineattrs = (color = red pattern = solid thickness = 0.5)
markerattrs=(color = red symbol=circlefilled);
run;
/* 模型初步识别 */
proc arima data=air;
identify var=plog(1) nlag=30;/*plog(1)即是进行一阶差分，nlag表示滞后的阶数*/
run;
/* 模型定阶 */
proc arima data=air;
identify var=plog(1) nlag=30 minic p=(0:14) q=(0:14);
/* 还可以添加选项minic, esacf, scan */
run;
/* 模型拟合 */
proc arima data=air;
identify var=plog(1) nlag=30;
estimate p=(8,12) q=(1,3) method=cls;
/* 还可以添加选项method=ML（极大似然）、ULS（非条件最小二乘法）、CLS（最小二乘法） */
run;
/* 模型预测 */
proc arima data=air plots(only)=forecast(forecast) /*plots=all*/;
identify var=plog(1) nlag=30;
estimate p=(8,12) q=(1,3) method=cls;
/*不要使用id属性，因为id增加会按照最小时间间隔增加（这里是秒，不是按月），
这样画出的预测图上的多步预测值几乎都在挤在一起。*/
forecast lead=12 out=out;/*lead=12表示向前预测12步，并输出到数据out*/
run;
/* proc print data=out; */
/* title 'Results(no id)'; */
/* run; */
/* 模型预测(获得带id的数据) */
proc arima data=air;
identify var=plog(1) nlag=30;
estimate p=(8,12) q=(1,3) method=cls;
forecast lead=12 id=Month out=outAir;
run;
/* proc print data=outAir; */
/* title 'Results(have id)'; */
/* run; */
/* 将实际数据和预测数据放在一张图上做对比 */
proc sgplot data=outAir;
series x=Month y=plog / legendlabel = 'actual'
markers lineattrs = (color=black pattern=dash thickness = 0.5)
markerattrs=(color=black symbol=diamond);
series x=Month y=forecast / legendlabel = 'forecast'
markers lineattrs = (color=red pattern=solid thickness = 0.5)
markerattrs=(color=red symbol=dot);
/* 置信区间为95%的下界 */
/* series x=Month y=l95 / legendlabel = 'l95' */
/* markers lineattrs = (color=green pattern=1 thickness = 0.5) */
/* markerattrs=(color=green symbol=star); */
/* 置信区间为95%的上界 */
/* series x=Month y=u95 / legendlabel = 'u95' */
/* markers lineattrs = (color=blue thickness = 0.5) */
/* markerattrs=(color=blue symbol=circlefield); */
title'Air Passengers';
run;
（2）AirPassengers.csv内容
Month,#Passengers
1949-01,112
1949-02,118
1949-03,132
1949-04,129
1949-05,121
1949-06,135
1949-07,148
1949-08,148
1949-09,136
1949-10,119
1949-11,104
1949-12,118
1950-01,115
1950-02,126
1950-03,141
1950-04,135
1950-05,125
1950-06,149
1950-07,170
1950-08,170
1950-09,158
1950-10,133
1950-11,114
1950-12,140
1951-01,145
1951-02,150
1951-03,178
1951-04,163
1951-05,172
1951-06,178
1951-07,199
1951-08,199
1951-09,184
1951-10,162
1951-11,146
1951-12,166
1952-01,171
1952-02,180
1952-03,193
1952-04,181
1952-05,183
1952-06,218
1952-07,230
1952-08,242
1952-09,209
1952-10,191
1952-11,172
1952-12,194
1953-01,196
1953-02,196
1953-03,236
1953-04,235
1953-05,229
1953-06,243
1953-07,264
1953-08,272
1953-09,237
1953-10,211
1953-11,180
1953-12,201
1954-01,204
1954-02,188
1954-03,235
1954-04,227
1954-05,234
1954-06,264
1954-07,302
1954-08,293
1954-09,259
1954-10,229
1954-11,203
1954-12,229
1955-01,242
1955-02,233
1955-03,267
1955-04,269
1955-05,270
1955-06,315
1955-07,364
1955-08,347
1955-09,312
1955-10,274
1955-11,237
1955-12,278
1956-01,284
1956-02,277
1956-03,317
1956-04,313
1956-05,318
1956-06,374
1956-07,413
1956-08,405
1956-09,355
1956-10,306
1956-11,271
1956-12,306
1957-01,315
1957-02,301
1957-03,356
1957-04,348
1957-05,355
1957-06,422
1957-07,465
1957-08,467
1957-09,404
1957-10,347
1957-11,305
1957-12,336
1958-01,340
1958-02,318
1958-03,362
1958-04,348
1958-05,363
1958-06,435
1958-07,491
1958-08,505
1958-09,404
1958-10,359
1958-11,310
1958-12,337
1959-01,360
1959-02,342
1959-03,406
1959-04,396
1959-05,420
1959-06,472
1959-07,548
1959-08,559
1959-09,463
1959-10,407
1959-11,362
1959-12,405
1960-01,417
1960-02,391
1960-03,419
1960-04,461
1960-05,472
1960-06,535
1960-07,622
1960-08,606
1960-09,508
1960-10,461
1960-11,390
1960-12,432
