[转载]教你用MATLAB做季节性调整

原文地址：教你用MATLAB做季节性调整作者：盘后TeaTime

时间序列分解的含义是，将时间序列分解成若干个因子叠加的共同影响。一般认为时间序列可以分解成三个因素：趋势-循环因素、季节性因素和白噪音。按照组合方式的不同，可以分为加法模型、乘法模型和对数加法模型。通常把去除季节性因素的过程叫做季节性调整。

 

季节调整的意义

 

由于生产活动中一般都要提前安排来年或者下个月甚至下周的生产采购和销售计划。所以对当前整体经济活动的趋势判断就显得至关重要。对于年度序列的数据而言，由于反映的是一个完整年份的经济运行状况，因此只要比较历史数据一般就可以得到每年的趋势变化情况。可是要得到一整年的数据情况才能判断趋势的话，却往往会出现落差很大或者超出预期的情况。比如提前观察到下个月经济开始收缩，但是企业反而根据去年制订的计划加快生产，显然就不合时宜。

 

正所谓计划没有变化快，所以就有了每个季度或者每月的公布数据，有助于提前进行对计划的修正。但是使用季度或者月度数据就会遇到季节效应的问题。比如说夏季的时候空调销量很好，可是如果不能剔除掉季节性影响的话，单纯与上月相比，很容易判断成销售加快，趋势上升。所以为了克服季节性因素对数据的影响，因此就必须使用季节性调整办法来对数据进行修正。当然更理想化来说，我们希望得到趋势因素。

 

国家统计局公布的数据一般是同比的。虽然说与去年同期比较，相对于未经季调的环比数据有优势。但是它不能及时判断出拐点。有研究表明，使用未经处理的同比数据对经济周期判断往往滞后于实际的转折点。而且我们国家还有个“特色”，很多数据都是1、2月合计公布的，统计局说这是为了避免1月份和2月份数据产生较大误差的处理方式。但实际上嘛，个人认为当所有人都明白1月和2月数据是受到季节性影响的时候，市场就不会对此有太大反应。反而容易让人觉得，统计局找这么个理由完全是为了省点工作。当然我没法改变统计局公布数据的方式和计划，但是实在要吐吐苦水，不公布单独的一月和二月数据对我们这些要做分析数据的研究人员实在有点为难。

 

关于Matlab的季节性调整方法

 

唠叨那么多，回归到正题。Matlab没有专门做季节调整的函数，这是让我觉得不解的地方。也许没有人会有这个需求，因为像Eviews这些软件都有相对应的工具，按个键就好了，非常方便。但是在使用这些季节性调整工具软件的时候，很容易就会犯错误。因为对各个参数设定不了解，同时也不能按照中国的情况来进行模型修改，所以直接用来做季调跟自己说不清楚，也没法跟别人解释清楚。中国人民银行有自己开发的一套软件，叫做中国版的X-12-ARIMA，也出了一套书（文献[a]）去说明这个软件。当然这些都要花大价钱买，对于广大的屌丝研究员而言，根本没有机会去用。

 

在Masthworks的网站上可以找到季节性调整相关的概念，以及季节性调整的几个例子。其中比较详尽的是利用复合移动平均法做季调的例子。[b] 这个例子基本涵盖了其余与季调相关的过程。如果对季调有所了解的朋友，应该能看出来这个例子是不考虑对极端情况以及节假日效应处理的估计，单纯连续使用移动平均法对数据进行季调，也就是经常用到的X-11方法。

 

关于X-11做季调的流程步骤，可以参考下表（参考[1] P134）

 

表1 X-11计算原型（乘法模型）

---

使用MATLAB对CPI环比进行季节性调整

 

准备数据

 

本文选用未经季调过的CPI环比数据来做展示。这个数据开始于1995年的1月份，到2013年7月结束。把1994年12月31日定义为基期，基期指数设为100。将所有环比数据还原为CPI的指数序列，得到图1。

 

图1 CPI走势图

 

从图上看，CPI有明显的上升趋势。但是季节性因素却不明显，仅有在1997年至2003年这段时间才呈现明显的季节性规律，主要是由于上升趋势较为明显。

 

1、使用2×12移动平均估计趋势-循环成分

 

由上图来看，使用乘法模型对该数据建模较合适。为了去除趋势成分，参考Matlab给出的例子，使用S2×13的中心对称移动平均对其进行平滑处理以估计趋势-循环的部分。但是中心化的移动平均容易丢失原数据序列两端的数据信息。

 

这段程序清楚展示了如何利用S2×13中心对称移动平均进行平滑。由于平滑处理后的数据首尾两端的信息会被丢失，此处的简单处理办法是使首尾6个数据等于邻近的数据值。

%% 1、使用2×12移动平均估计趋势-循环成分：

sW13 = [1/24;repmat(1/12,11,1);1/24];
CPIs = conv(CPI,sW13,'same');
CPIs(1:6) = CPIs(7);
CPIs(N-5:N) = CPIs(N-6);

可以将平滑处理后的数据与原始数据作比较，在图2上用红线表示的是平滑处理后的数据。可以看到首尾是水平无变化的。

 

图2 初步平滑后的序列与原始序列对比

 

2、估计季节-不规则成分：

 

用原始数据除以平滑处理后的数据就得到去除趋势成分的CPI数据，它包含了周期成分和不规则的成分（白噪音）。在这个基础上可以进一步分离出周期成分，也就是季节性的成分。这只是第一步，并非是最后得到的季节成分。后面还要反复运用移动平均来做调整。

 

 

%% 2、估计季节-不规则成分：
xt = CPI./CPIs;

 

3、对每个月份应用3×3移动平均估计预备季节成分

 

Matlab要输入S3×3阶的移动平均算子必须手动输入。

% S3x3 季调滤子
% 中心对称系数
sW3 = [1/9;2/9;1/3;2/9;1/9];
% 非中心对称系数（用于调整首尾6个数据）
aW3 = [.259 .407;.37 .407;.259 .185;.111 0];

aW3是Musgrave（1964a,1964b）给出的用于估计季节因子的复合移动平均相对应的非对称滤子，但是这在X-11方法中并未被采用。对于这些滤子的合理性，尚未看到有任何出版物对其作出解释。[a] 但是Matlab里面却调用了这些滤子，在实际运用中需要当心。

% 使用2×12移动平均估计预备季节成分
sW13 = [1/24;repmat(1/12,11,1);1/24];
sb = conv(shat,sW13,'same');
sb(1:6) = sb(s+1:s+6);
sb(N-5:N) = sb(N-s-5:N-s);

% 然后标准化
s33 = shat./sb;

 

图3 预备季调因子

 

4、估计季节调整后序列

得到预备季调因子以后，便可得到第一阶段季调后的序列

dt = CPI./s33;

 

5、用13项Henderson移动平均估计趋势-循环成分

在X-11中，使用Henderson移动平均从季调后的序列中提取趋势-循环成分。本例中使用的是13项Henderson移动平均算子

 

% Henderson 系数
sWH = [-0.019;-0.028;0;.066;.147;.214;
      .24;.214;.147;.066;0;-0.028;-0.019];

 

对应的非对称移动平均系数。但对比[a]中所给的系数略有不同。

 

% 非对称系数（针对首位6个数据）
aWH = [-.034  -.017   .045   .148   .279   .421;
       -.005   .051   .130   .215   .292   .353;
        .061   .135   .201   .241   .254   .244;
        .144   .205   .230   .216   .174   .120;
        .211   .233   .208   .149   .080   .012;
        .238   .210   .144   .068   .002  -.058;
        .213   .146   .066   .003  -.039  -.092;
        .147   .066   .004  -.025  -.042  0    ;
        .066   .003  -.020  -.016  0      0    ;
        .001  -.022  -.008  0      0      0    ;
       -.026  -.011   0     0      0      0    ;
       -.016   0      0     0      0      0    ];

 

应用这些滤子便可以得到最终的趋势成分h13

 

% h13就是最终的趋势成分
first = 1:12;
last = N-11:N;
h13 = conv(dt,sWH,'same');
h13(N-5:end) = conv2(dt(last),1,aWH,'valid');
h13(1:6) = conv2(dt(first),1,rot90(aWH,2),'valid');

 

对应的季节成分

 

 

% 估计季节-不规则成分：
xt = CPI./h13;

 

图4 13项Henderson移动平均法得到的趋势因素

 

6、对每个月份应用3×5移动平均估计最终季节成分

 

% S3x5 季调滤子
% 中心对称系数
sW5 = [1/15;2/15;repmat(1/5,3,1);2/15;1/15];
% 非中心对称系数
aW5 = [.150 .250 .293;
       .217 .250 .283;
       .217 .250 .283;
       .217 .183 .150;
       .133 .067    0;
       .067   0     0]; 

……

% 运用2×12移动平均
sW13 = [1/24;repmat(1/12,11,1);1/24];
sb = conv(shat,sW13,'same');
sb(1:6) = sb(s+1:s+6);
sb(N-5:N) = sb(N-s-5:N-s);

% 然后标准化得到最终季节成分
s35 = shat./sb;
Season = s35;

 

7、估计季节调整后序列

用原始数据除以季节因子就可以得到最终季调后的序列

 

 

 dt = CPI./s35;

图5 季调后CPI走势

 

8、估计最终趋势

最终趋势便是前面的h13，这里给出各个成分序列的对比图。

 

图6 CPI各成分序列的对比图

 

9、估计最终不规则成分

用季调后的序列除以趋势序列便得到不规则序列

 

 

Irr = dt./h13;

 

图7 CPI的不规则成分（白噪音）

 

---

后记

 

1、本文主要参考资料有二：

  a)中国人民银行调查统计司编著的《时间序列X-12-ARIMA季节调整——原理与方法》，其中对X-12-ARIMA做季调的程序做了非常详细的说明。至今统计局仍采用这一程序进行季节性调整。

  b)MATHWORKS主页上关于季调的文档。例子中做季调的对象是航空公司乘客统计的月度数据。此处的代码都基本出自于这个例子的源程序。

http://www.mathworks.cn/cn/help/econ/seasonal-adjustment-using-snxd7m-seasonal-filters.html

 

2、本文使用的X-11方法一直饱受诟病，对它的批评集中在下面三点：缺乏统计理论支撑，有时候会扭曲各种构成成分，会扭曲变量之间的关系。后来Box和Jenkins（1970）有关ARIMA模型的成果对基于模型的季节调整方法做出了重大贡献，后期形成了X-11-ARIMA方法。后来再在前者基础上出现了X-12-ARIMA方法，基本思路都是在使用X-11进行季调前必须对序列进行建模，以弥补X-11方法的缺陷。

 

3、往后会介绍ARIMA模型，以及如何结合ARIMA进行季调的方法。当然未来可能还会涉及到异常值处理以及节假日调整等问题。这个系统比较庞大，所以本文并不打算就此展开。后期将逐一介绍，请关注。