频域线性回归

打开实时脚本

这个例子展示了如何使用离散傅里叶变换来构造时间序列的线性回归模型。本例中使用的时间序列是1973年至1979年美国每月意外死亡人数。数据发表在Brockwell and Davis(2006)上。原始来源是美国国家安全委员会。

输入数据。复制exdata矩阵进入MATLAB®工作空间。

Exdata = [9007 7750 8162 7717 7792 7836 8106 6981 7306 7461 6957 6892 8928 8038 8124 7776 7870 7925 8106 8129 10017 9387 8634 8890 9115 10120 10093 10078 10625 10484 10744 9823 9620 9179 9302 9827 9713 8743 8285 8037 8314 9110 9938 9129 8433 8488 8850 9070 9161 8710 8160 7874 8647 8796 9240];

exdata是一个12 × 6矩阵。的每一列exdata包含12个月的数据。每一列的第一行包含相应年份1月份美国意外死亡人数。每一列的最后一行包含相应年份12月美国意外死亡人数。

将数据矩阵重塑为72 × 1的时间序列，并绘制出1973年至1978年的数据。

Ts =重塑(exdata,72,1);年= linspace(1973,1979,72);情节(年,ts,“啊——”，“MarkerFaceColor”，“汽车”)包含(“年”) ylabel (意外死亡人数）

图中包含一个轴对象。axis对象包含一个line类型的对象。

对数据的目视检查表明，意外死亡人数呈周期性变化。振荡的周期大约为1年(12个月)。数据的周期性表明，一个合适的模型可能是

$X （ n ）＝ μ + \sum_{k} （ {一个}_{k} 因为 \frac{2 π k n}{N} + B_{k} 因为 \frac{2 π k n}{N} ） + ε （ n ），$

在哪里 $μ$ 是总体平均值， $N$ 是时间序列的长度，和 $ε （ n ）$ 是均值为零且有一定方差的独立等分布(IID)高斯随机变量的白噪声序列。附加噪声项解释了数据中固有的随机性。模型的参数是总体均值和余弦和正弦的幅值。该模型的参数是线性的。

要在时域中构建线性回归模型，必须指定余弦和正弦的频率，形成设计矩阵，并求解正态方程，以获得模型参数的最小二乘估计。在这种情况下，更容易使用离散傅里叶变换来检测周期性，只保留傅里叶系数的子集，并反转变换以获得拟合的时间序列。

对数据进行频谱分析，以揭示哪些频率对数据的可变性有重要影响。由于信号的总体平均值约为9,000，并与0频率的傅里叶变换成正比，因此在频谱分析之前减去平均值。这降低了0频率的大幅度傅里叶系数，使任何显著的振荡更容易被检测到。傅里叶变换中的频率间隔为时间序列长度的倒数，1/72。数据按月采样，频谱分析的最高频率为1个周期/2个月。在这种情况下，可以方便地以周期/年为单位查看频谱分析，因此可以相应地缩放频率以进行可视化。

TSDFT = fft(ts-mean(ts));频率= 0:1/72:1/2;情节(频率。* 12、abs (tsdft(1:长度(ts) / 2 + 1)),“啊——”，.．.“MarkerFaceColor”，“汽车”)包含(“周期/年”) ylabel (“级”) ax = gca;斧子。XTick = [1/6 1 2 3 4 5 6];

图中包含一个轴对象。axis对象包含一个line类型的对象。

根据震级，1个周期/12个月的频率是数据中最显著的振荡。1个周期/12个月的幅度是其他幅度的两倍多。然而，光谱分析表明，数据中还存在其他周期性成分。例如，在1个周期/12个月的谐波(整数倍)上出现周期性成分。似乎还有一个周期成分，周期为1个周期/72个月。

根据数据的频谱分析，使用频率为最显著成分的余弦和正弦项拟合一个简单的线性回归模型:1个周期/年(1个周期/12个月)。

确定离散傅里叶变换中对应于1个周期/12个月的频率仓。因为频率间隔为1/72，第一个箱子对应于0频率，正确的箱子是72/12+1。这是频率箱积极的频率。还必须包括对应的频率仓负频率:-1周期/12个月。通过MATLAB索引，得到负频的频仓为72-72/12+1。

创建一个72乘1的零向量。用对应于1个周期/12个月的正频率和负频率的傅里叶系数填充向量的适当元素。对傅里叶变换进行逆变换，并将总体均值相加，以获得对意外死亡数据的拟合。

Freqbin = 72/12;freqbin = [freqbin 72-freqbin]+1;tfit = 0 (72,1);Tsfit (freqbins) = tsdft(freqbins);tfit = ifft(tfit);Mu = mean(ts);tfit = mu+ tfit;

用两个傅立叶系数将原始数据与拟合级数一起绘制出来。

情节(年,ts,“啊——”，“MarkerFaceColor”，“汽车”)包含(“年”) ylabel (意外死亡人数)举行在tsfit情节(几年,“线宽”2)传说(“数据”，“安装模式”)举行从

图中包含一个轴对象。axis对象包含2个line类型的对象。这些对象表示数据，拟合模型。

拟合模型似乎捕捉了数据的一般周期性，并支持数据以1年为周期振荡的初步结论。金宝app

为了评估1个周期/12个月的单一频率如何充分地解释观测到的时间序列，形成残差。如果残差类似于白噪声序列，则具有一个频率的简单线性模型已经充分模拟了时间序列。

为了评估残差，使用白噪声95%置信区间的自相关序列。

残留= ts- tfit;[xc，滞后]= xcorr(残余，50，多项式系数的）;茎(滞后(51:结束),我(51:结束),“填充”)举行在Lconf = -1.96*ones(51,1)/√(72);Uconf = 1.96*ones(51,1)/√(72);情节(lconf滞后(51:结束),“r”)情节(滞后(51:结束),uconf,“r”)包含(“滞后”) ylabel (“相关系数”)标题(“残差的自相关”)举行从

图中包含一个轴对象。标题为残差自相关的轴对象包含茎型、行型3个对象。

自相关值在一些滞后时落在95%置信范围之外。残差似乎不是白噪声。结论是，具有一个正弦分量的简单线性模型不能解释意外死亡数量中的所有振荡。这是意料之中的，因为光谱分析揭示了除了主导振荡之外的其他周期成分。创建一个包含光谱分析所指出的额外周期项的模型将改善拟合并使残差变白。

拟合一个由三个最大的傅里叶系数组成的模型。因为你必须保留对应于负频率和正频率的傅里叶系数，保留最大的6个指标。

Tsfit2dft = 0 (72,1);[Y,I] = sort(abs(tsdft)，“下”）;指数= I(1:6);Tsfit2dft(指数)= tsdft(指数);

证明仅保留72个傅里叶系数中的6个(3个频率)就能保留大部分信号的能量。首先，证明保留所有的傅里叶系数在原始信号和傅里叶变换之间产生能量等价。

规范(1 /√(72)* tsdft, 2) /规范(ts-mean (ts), 2)

Ans = 1.0000

比值是1。现在，检查只保留3个频率的能量比。

规范(1 /√(72)* tsfit2dft, 2) /规范(ts-mean (ts), 2)

Ans = 0.8991

几乎90%的能量都被保留了下来。等效地，时间序列的90%的方差由3个频率分量解释。

根据3个频率分量对数据进行估计。将原始数据、一频率模型和三频率模型进行比较。

Tsfit2 = mu+ifft(tsfit2dft，“对称”）;情节(年,ts,“啊——”，“markerfacecolor”，“汽车”)包含(“年”) ylabel (意外死亡人数)举行在tsfit情节(几年,“线宽”2)图(年,tsfit2,“线宽”2)传说(“数据”，“1频率”，“3频率”)举行从

图中包含一个轴对象。axis对象包含3个line类型的对象。这些对象表示数据，1频率，3频率。

使用3个频率提高了对原始信号的拟合。你可以通过检查三频率模型残差的自相关来看到这一点。

Resid = ts-tsfit2;[xc，滞后]= xcorr(残余，50，多项式系数的）;茎(滞后(51:结束),我(51:结束),“填充”)举行在Lconf = -1.96*ones(51,1)/√(72);Uconf = 1.96*ones(51,1)/√(72);情节(lconf滞后(51:结束),“r”)情节(滞后(51:结束),uconf,“r”)包含(“滞后”) ylabel (“相关系数”)标题(“残差的自相关”)举行从

图中包含一个轴对象。标题为残差自相关的轴对象包含茎型、行型3个对象。

使用3个频率导致残差更接近白噪声过程。

证明从傅里叶变换得到的参数值等价于时域线性回归模型。通过形成设计矩阵并求解正态方程，找到三个频率的总体平均值、余弦振幅和正弦振幅的最小二乘估计。将拟合的时间序列与傅里叶变换得到的时间序列进行比较。

X = ones(72,7);X(:，2) = cos(2*pi/72*(0:71))';X(:，3) = sin(2*pi/72*(0:71))';X(:，4) = cos(2*pi*6/72*(0:71))';X(:，5) = sin(2*pi*6/72*(0:71))';X(:，6) = cos(2*pi*12/72*(0:71))';X(:，7) = sin(2*pi*12/72*(0:71))';beta = X\ts;tsfit_lm = X*beta;马克斯(abs (tsfit_lm-tsfit2))

Ans = 9.0949e-12

这两种方法得到的结果是相同的。两种波形之间差异的最大绝对值在10-12的数量级上。在这种情况下，频域方法比等效的时域方法更容易。您自然会使用频谱分析来直观地检查数据中存在哪些振荡。从这一步开始，使用傅里叶系数来构造一个由余弦和正弦组成的信号模型是很简单的。

有关时间序列频谱分析和时域回归的等效性的更多细节，请参见(Shumway和Stoffer, 2006)。

虽然光谱分析可以回答哪些周期性成分对数据的可变性有显著贡献，但它不能解释为什么这些成分会出现。如果仔细检查这些数据，就会发现12个月周期中的最小值往往出现在2月，而最大值则出现在7月。对这些数据的一个合理解释是，人们在夏天比冬天更活跃。不幸的是，由于这种活动的增加，发生致命事故的可能性也增加了。

参考文献

布洛克威尔，彼得·J·和理查德·a·戴维斯。时间序列:理论与方法．纽约:施普林格，2006。

沙姆韦，罗伯特H.和大卫S.斯托弗。时间序列分析及其应用(附R例)．纽约:施普林格，2006。

另请参阅

fft|传输线|xcorr