自协方差函数的Matlab实现

图灵猫

已于 2022-03-22 15:39:42 修改

阅读量3.7k

点赞数 1

文章标签：时间序列 matlab

于 2022-03-20 19:58:53 首次发布

本文链接：https://blog.csdn.net/qq_42818403/article/details/123619891

版权

自协方差函数的Matlab实现

文章目录

自协方差函数的Matlab实现

自协方差函数简介

对于满足均值遍历性和二阶矩遍历性的平稳时间序列一次具体观测值 $\left\{x_{t}, t=1,2, \cdots, N\right\}$ , 总体平均可转化为时间平均. 因此, 样本自协方差函数也通常可由下面两种方式计算得到:
$\begin{cases}\hat{\gamma}_{k}=\frac{1}{N} \sum_{t=1}^{N-k}\left(x_{t}-\hat{\mu}\right)\left(x_{t+k}-\hat{\mu}\right), & 0 \leqslant k \leqslant N-1 \\ \hat{\gamma}_{k}=\hat{\gamma}_{-k}, & 1-N \leqslant k \leqslant-1,\end{cases}\tag{1}$
或
$\begin{cases}\hat{\gamma}_{k}=\frac{1}{N-k} \sum_{t=1}^{N-k}\left(x_{t}-\hat{\mu}\right)\left(x_{t+k}-\hat{\mu}\right), & 0 \leqslant k \leqslant N-1 \\ \hat{\gamma}_{k}=\hat{\gamma}_{-k}, & 1-N \leqslant k \leqslant-1\end{cases}\tag{2}$
其中, $\hat{\mu}$ 为平稳序列的均值的估计值. 在实际应用中, 样本自协方差函数的计算通常采用第一种形式.

其原因在于, 利用 $(1)$ 式定义的样本自协方差函数能使得 $N$ 阶的样本自协方差矩阵
$\hat{\boldsymbol{\Gamma}}_{N}=\left(\hat{\gamma}_{k-j}\right), k, j=1,2, \cdots, N$
正定.

为了计算自协方差矩阵(ACM),一种最直接的方式是利用两层循环嵌套结构计算,但此方式过于“原始”,对于长期序列计算速度相对比较慢.

注意到, 只要
$y_{j}=x_{j}-\hat{\mu}, \quad j=1,2, \cdots, N$
不全为零, 矩阵
$\boldsymbol{A}=\left(\begin{array}{cccccccc} 0 & \cdots & 0 & y_{1} & y_{2} & \cdots & y_{N-1} & y_{N} \\ 0 & \cdots & y_{1} & y_{2} & y_{3} & \cdots & y_{N} & 0 \\ \vdots & & \vdots & \vdots & \vdots & & \vdots & \vdots \\ y_{1} & \cdots & y_{N-1} & y_{N} & 0 & \cdots & 0 & 0 \end{array}\right)_{N \times(2 N-1)}$
是行满秩的. 样本自协方差矩阵可表示为
$\hat{\boldsymbol{\Gamma}}_{N}=\left(\begin{array}{cccc} \hat{\gamma}_{0} & \hat{\gamma}_{1} & \cdots & \hat{\gamma}_{N-1} \\ \hat{\gamma}_{1} & \hat{\gamma}_{0} & \cdots & \hat{\gamma}_{N-2} \\ \vdots & \vdots & & \vdots \\ \hat{\gamma}_{N-1} & \hat{\gamma}_{N-2} & \cdots & \hat{\gamma}_{0} \end{array}\right)=\frac{1}{2} \boldsymbol{A} \boldsymbol{A}^{\mathrm{T}}$

于是,得益于计算机对矩阵并行运算的优化.我们可以采用数组化编程思想,利用矩阵运算,轻松且快速的计算出自协方差矩阵.

上机实验报告

数据集

为了感受不同时间序列的统计性质,选用以下两种数据集分别建模.

洗发水销售数据集:shampoo.csv
上证指数数据集:SH600031.csv

实验过程

洗发水销售数据集

该数据集描述了每月洗发水的销售量

单位是销售计数，有36个观察值.

分布图

程序源代码

filename='shampoo.csv';
data = readmatrix(filename, 'OutputType', 'single');
x=data(1:36,2);
y=x-mean(x);
A=zeros(length(y),2*length(y)-1);
for i=1:length(y)
    A(i,length(y)-i+1:length(y)-i+length(y))=y;
end
ACM=1/length(y).*A*A';

运行结果分析:
每月销售量时序图

自协方差函数

自协方差函数最大值:2.156610448813604e+04

自协方差函数最小值:-6.896521602231956e+03

从图像及数据可以发现:洗发水月销售量具有较明显的子相关性.

股票信息数据集

该数据集描述了上证指数每日的开盘，最高，最低，收盘，成交量，成交额.
取连续365天的开盘价进行建模.

分布图

程序源代码

filename='SH600031.csv';
data = readmatrix(filename, 'OutputType', 'single');
x=data(1:365,3);
y=x-mean(x);
A=zeros(length(y),2*length(y)-1);
for i=1:length(y)
    A(i,length(y)-i+1:length(y)-i+length(y))=y;
end
ACM=1/length(y).*A*A';

运行结果分析:
每日上证指数开盘价时序图