基于Box-Cox变换的改进ARMA模型的MATLAB实现

最新推荐文章于 2024-05-14 21:49:33 发布

三只佩奇不结义

最新推荐文章于 2024-05-14 21:49:33 发布

阅读量3.2k

点赞数 8

分类专栏： MATLAB随笔文章标签：算法 matlab

本文链接：https://blog.csdn.net/qq_41196612/article/details/104962315

版权

MATLAB随笔专栏收录该内容

15 篇文章 10 订阅

订阅专栏

文章目录

1 问题背景
2 Box-Cox变换
- 2.1 介绍
- 2.2 极大似然参数估计
3 ARMA模型的建立
4 代码

前言： 第一次听说Box-Cox变换，发现原来这玩意这么厉害，值得写一写记录一下。

以下内容参考了栗然等人写的论文《基于Box-Cox变换的风电场短期风速预测模型》

1 问题背景

风力发电是解决当前突出的能源和环境问题的有效手段，因而得到了普遍重视和大规模开发利用，是目前世界上增长速度最快的能源开发形式。尤其近十几年来其发电成本大幅下降，已具备与传统常规电源竞争的潜力。

目前，国内外对于风力发电的研究越来越深入和广泛，但关于风电场风速预测的研究还不能达到令人满意的程度。风电场风速预测的方法较多，最简单的预测方法是持续法，即把最近一点的风速观测值作为下一点的预测值；其他预测方法有卡尔曼滤波法(Kalman filters)、时间序列法(ARMA)、人工神经网络法(ANN)、模糊逻辑法(Fuzzy Logic)、空间相关法等。时间序列法作为短期风速预测的方法之一，只需知道风电场的单一风速时间序列即可建立预测模型，且预测效果较好。但传统的时间序列ARMA $(p, q)$ 预测模型需要寻找 $y_t$ 与 $y_{t-i}$ 、 $a_{t-j}$ 的线性关系，即：

$y_{t}=\sum_{i=1}^{p} \varphi_{i} y_{t-i}+a_{t}-\sum_{j=1}^{q} \theta_{j} a_{t-j}$

式中， $\varphi_{i}(i=l, 2, \cdots, p)$ 为自回归参数； $\theta_{j}(j=1,2,\cdots,q)$ 为滑动平均参数； $a_t$ 表示随机干扰，是一零均值、方差为 $\sigma_{a}^{2}$ 的正态白噪声过程，即 $a_t \sim N(0,\sigma_{a}^2)$ 而实际的风速数据并不总是满足以上关系，Box-Cox变换能够在一定程度上解决该问题。它通过对因变量的变换，使得变换后的因变量与自变量具有线性关系，误差各分量是等方差且相互独立的。

本文以某风电场为例，首先通过Box-Cox变换对原始风速时间序列进行转换，使用极大似然法，确定风速目标值的最优转换形式，满足ARMA $(p, q)$ 模型建立条件；针对转换后的风速数据，利用BIC准则对ARMA(P，q)模型定阶，确定合适阶数模型进行风电场风速预测。结果分析表明，该方法误差较小，预测曲线能很好地跟踪风速变化趋势，有效提高了风速预测精度。

2 Box-Cox变换

2.1 介绍

在做线性回归的过程中，一般线性模型假定：

$\beta+\varepsilon, \varepsilon \sim N\left(0, \delta^{2} I\right)$

该模型需要满足四个条件：

线性性： $E (Y)$ 是 $X$ 中各变量的线性函数；
独立性： $\varepsilon_1,\varepsilon_2,\cdots,\varepsilon_n$ 相互独立；
方差齐性： $D(\varepsilon_1)=D(\varepsilon_2)=\cdots=D(\varepsilon_n)=\delta^2$
正态性： $\varepsilon_1,\varepsilon_2,\cdots,\varepsilon_n$ 服从正态分布

在处理实际经济问题和社会问题时，由于海量数据比较凌乱，同时在建立回归模型时，个别变量的系数通不过。例如生物医学等数据的特殊性，往往不可观测的误差 $\varepsilon$ 可能是和预测变量相关的，不服从正态分布，于是给线性回归的最小二乘估计系数的结果带来误差，为了满足上述四个条件而不丢失信息，有时需要改变一下数据形式，进而Box-Cox变换得到了广泛推广。

Box和Cox于1964年提出了一种基于极大似然法的幂转换模型。Box-Cox幂分布族是一种十分有用的连续分布族。其转换模型为

$y^{\lambda}=\left\{\begin{array}{ll} \frac{y^{\lambda}-1}{\lambda}, & \lambda \neq 0 \\ \ln y, & \lambda=0 \end{array}\right.$
这里 $\lambda$ 是一个待定的变换参数。对不同的 $\lambda$ ，所做的变换自然不同，所以就是一个变换族，当 $\lambda=0$ 时，该变换为对数变换。该变换的逆变换为
$y=\left\{\begin{array}{l} \left(1+\lambda y^{(\lambda)}\right)^{\frac{1}{\lambda}}, \lambda \neq 0 \\ \exp \left(y^{(\lambda)}\right), \lambda=0 \end{array}\right.$

但是在变换前可能用到 $l o g (y)$ ，因此必须要求 $y > 0$ 。对于一般的数据，对于任意取值的 $y$ 的Box-Cox变换可用下面的公式
$y^{(\lambda)}=\left\{\begin{array}{l} \frac{(y+c)^{\lambda}-1}{\lambda g}, \lambda \neq 0 \\ \frac{\log (y+c)}{g}, \lambda=0 \end{array}\right.$

其中，参数 $c$ 就是为了使 $y + c > 0$ ， $g$ 的默认取值为1，也可取 $\hat{y}^{\lambda+1}$ ， $\hat{y}$ 为 $y$ 的几何平均。显然参数 $c$ 的取值很好确定，公式中需要估计的参数就是 $\lambda$ 。为了简便处理，一般理论上我们假定反应变量 $y > 0$ 。

有一个很典型的直观例子：
变换前
在这里插入图片描述
变换后

2.2 极大似然参数估计

对因变量的观察值 $y_1,\cdots,y_n$ ，应用上述变换，得到变换后的变量为

$y^{(\lambda)}=\left(y_{1}^{(\lambda)}, \cdots, y_{n}^{(\lambda)}\right)$
这就是说，要求通过因变量的变换，使得变换后的 $y^{(\lambda)}$ 与自变量具有线性依托关系。因此，Box-Cox变换是通过参数的适当选择，达到对原来数据的“综合治理”，使其满足一个线性模型条件。

对于 $\lambda$ 值的选择，可以通过极大似然法来估计。首先，在一个经验范围内选择参数 $\lambda$ 的值，然后使用下式计算。
$L^{(\lambda)}=-\frac{n}{2} \ln \sigma_{a}^{2}+\ln J(\lambda, y)$

式中，对于所有的 $\lambda$ ，有
$\ln J(\lambda, y)=\prod_{i=1}^{n} \frac{\partial W_{i}}{\partial y_{i}}=\prod_{i=1}^{n} y_{i}^{\lambda-1}$

对于每一个 $\lambda$ 来说，\sigma_{a}^{2}是 $y^{(\lambda)}$ 的极大似然估计，可通过下式计算得到。
$\sigma_{a}^{2}=\frac{1}{n} \sum_{i=1}^{n}\left(y_{i}^{(\lambda)}-\bar{y}^{(\lambda)}\right)^{2}$

经推导可得到如下方程：
$L^{(\lambda)}=-\frac{n}{2} \cdot \ln \left[\sum_{i=1}^{n} \frac{\left(y_{i}^{(\lambda)}-\bar{y}^{(\lambda)}\right)^{2}}{n}\right]+(\lambda-1) \cdot \sum_{i=1}^{n} \ln y_i$

式中， $\bar{y}^ {(\lambda)} = \frac{1}{n}\sum_{i=1}^{n}y_{i}^{(\lambda)}$

每一个 $\lambda$ 对应 $y^{(\lambda)}$ 的都可得到相应的 $L^{(\lambda)}$ 。由此可以描绘相应的 $\lambda$ 与 $L^{(\lambda)}$ 的关系图，从中我们可以得到相应的最优 $\lambda^{*}$ ，使得 $L^{(\lambda)}$ 最大；该优化的 $\lambda$ 对应了最优的转换模型。如下图风速序列Box-Cox转换分析曲线( $\lambda$ ∈[-2,2]，步长为0.01)所示，可以容易地找出 $\lambda^{*}$ 。在本文中找出的最优值为-0.2200。
在这里插入图片描述

3 ARMA模型的建立

ARMA模型不是本文的重点，所以就不多说了。在本文中使用了AIC确定ARMA $(p, q)$ 中的 $p$ 和 $q$ ，Box-Cox变换最终选取的模型为ARMA(3,3)，原始序列最终选取的模型为ARMA(2,3)。

利用建立的模型做预测，得到的结果如下图
在这里插入图片描述

4 代码

主程序

clc,clear

data = xlsread('新风速.xlsx');
y = data;
%% 极大似然法确定最优的lambda
lambda = -10:0.01:10;

L = zeros(1,length(lambda));
for i = 1:length(lambda)
    % 做Box-Cox变换
    if lambda(i) == 0
        y_lambda = log(y);
    else
        y_lambda = (y.^lambda(i)-1)/lambda(i);
    end
    
    % 计算似然函数
    y_hat = mean(y_lambda);
    n = length(y);
    L(i) = -n/2*log( ( sum((y_lambda - y_hat).^2) )/n ) + (lambda(i)-1)*sum(log(y));
end
plot(lambda,L)   % 似然函数随lambda的变化情况
xlabel('参数λ')
ylabel('似然函数L')

% 找到使得似然函数取得最大值的lambda
[~,ind] = max(L);
best_lambda = lambda(ind);

% 用最优的lambda最Box-Cox变换
if best_lambda == 0
    y_lambda = log(y);
else
    y_lambda = (y.^best_lambda-1)/best_lambda;
end

%% ARMA模型拟合与预测
n = 20;  % 预测未来数据的个数
% 经过Box-Cox变换的
[y_forecast1,best_RM1] = ARMA_function(y_lambda,n);

% 未经过变换的
[y_forecast2,best_RM2] = ARMA_function(y,n);
disp(['Box-Cox变换最终选取的模型为ARMA(' num2str(best_RM1(1)) ',' num2str(best_RM1(2)) ')'])
disp(['原始序列最终选取的模型为ARMA(' num2str(best_RM2(1)) ',' num2str(best_RM2(2)) ')'])

%% Box_Cox逆变换
if best_lambda == 0
    y_inv = exp(y_lambda);
    y_inv_fore = exp(y_forecast1);
else
    y_inv = (1 + best_lambda*y_lambda).^(1/best_lambda);
    y_inv_fore = (1 + best_lambda*y_forecast1).^(1/best_lambda);
end

%% 绘预测图
% 经过Box-Cox变换的
figure
m = length(y_inv);
y_inv = [y_inv;y_inv_fore(1)];
plot(1:m+1,y_inv)
hold on
plot(m+1:m+n,y_inv_fore)
title('经过Box-Cox处理后的预测情况')

% 未经过变换的
figure
m = length(y);
y = [y;y_forecast2(1)];
plot(1:m+1,y)
hold on
plot(m+1:m+n,y_forecast2)
title('原始序列预测情况')a

ARMA模型构建与预测的函数

function [y_forecast,best_RM] = ARMA_function(y_lambda,n)
%% 平稳性检验(pp检验看是否有单位根)
% 如果输出的hp为1，则说明该数据没有单位根是平稳数据，可以进行建模
disp('使用PP检验，如果不能拒绝原假设则说明指数序列存在单位根')
[hp,hpValue] = pptest(y_lambda,'model','TS');
hp

%% 模型定阶
k=0; %初始化试探模型的个数
for i=0:3
    for j=0:3
        if i==0 && j==0
            continue
        elseif i==0
            ToEstMd=arima('MALags',1:j,'Constant',0); %指定模型的结构
        elseif j==0
            ToEstMd=arima('ARLags',1:i,'Constant',0); %指定模型的结构
        else
            ToEstMd=arima('ARLags',1:i,'MALags',1:j,'Constant',0); %指定模型的结构
        end
        k=k+1; R(k)=i; M(k)=j;
        [~,EstParamCov,logL,info]=estimate(ToEstMd,y_lambda); %模型拟合
        numParams = sum(any(EstParamCov)); %计算拟合参数的个数
        %compute Akaike and Bayesian Information Criteria
        [aic(k),bic(k)]=aicbic(logL,numParams,length(y_lambda));
    end
end
% fprintf('R,M,AIC,BIC 的对应值如下\n '); %显示计算结果
% check = [R',M',aic',bic']

% 找到最小BIC对应的R和M值
[~,ind] = min(aic);
best_RM(1) = R(ind);
best_RM(2) = M(ind);
%% 模型预测
ToEstMd=arima('ARLags',1:best_RM(1),'MALags',1:best_RM(2),'Constant',0); %指定模型的结构
[EstMd,EstParamCov,logL,info] = estimate(ToEstMd,y_lambda); %模型拟合
y_forecast = forecast(EstMd,n,'Y0',y_lambda); %计算n步预报值,注意已知数据是列向量
end

本文使用的数据

10.0725000000000
9.33750000000000
9.77250000000000
9.55125000000000
10.1212500000000
9.36750000000000
9.27375000000000
10.3725000000000
9.92250000000000
10.0312500000000
10.1775000000000
9.05250000000000
9.22500000000000
10.0012500000000
9.18375000000000
9.52875000000000
9.87750000000000
9.48750000000000
10.3575000000000
9.68625000000000
10.2075000000000
10.1662500000000
9.94875000000000
9.89250000000000
9.41625000000000
10.3125000000000
9.85500000000000
10.1887500000000
10.1437500000000
9.03000000000000
9.67875000000000
9.64500000000000
9.95250000000000
10.0200000000000
9.48375000000000
9.90750000000000
9.86250000000000
9.81375000000000
10.3312500000000
10.9950000000000
9.71250000000000
10.4775000000000
9.67875000000000
10.3950000000000
10.8937500000000
10.2112500000000
10.6950000000000
9.99375000000000
10.9237500000000
11.6662500000000
9.50625000000000
10.9987500000000
11.7412500000000
11.8650000000000
11.8537500000000
11.4300000000000
11.6025000000000
11.5350000000000
10.9462500000000
10.2600000000000
9.79125000000000
10.3762500000000
10.0537500000000
9.98250000000000
11.2762500000000
10.0012500000000
10.5600000000000
10.4662500000000
11.1037500000000
11.0287500000000
11.4712500000000
11.1262500000000
10.6500000000000
10.8787500000000
11.2350000000000
10.6425000000000
10.2187500000000
11.6175000000000
11.4187500000000
10.5262500000000
10.9687500000000
10.5450000000000
10.0500000000000
11.1300000000000
11.5425000000000
12.1275000000000
11.5800000000000
11.2462500000000
11.4262500000000
12.2437500000000
11.8725000000000
11.2650000000000
12.2625000000000
11.8050000000000
10.8712500000000
11.0362500000000
11.5500000000000
11.0287500000000
10.8000000000000
11.1825000000000
10.0762500000000
11.0400000000000
10.5037500000000
10.9237500000000
10.8112500000000
10.4250000000000
11.2237500000000
10.4662500000000
11.3700000000000
10.1212500000000
10.7962500000000
10.9575000000000
10.5262500000000
10.6200000000000
10.8262500000000
11.0212500000000
10.5712500000000
10.2225000000000
10.8862500000000

三只佩奇不结义

关注

8
点赞
踩
67

收藏

觉得还不错? 一键收藏
打赏
2
评论
基于Box-Cox变换的改进ARMA模型的MATLAB实现

前言：第一次听说Box-Cox变换，发现原来这玩意这么厉害，值得写一写记录一下。以下内容参考了栗然等人写的论文《基于Box-Cox变换的风电场短期风速预测模型》1 问题背景风力发电是解决当前突出的能源和环境问题的有效手段，因而得到了普遍重视和大规模开发利用，是目前世界上增长速度最快的能源开发形式。尤其近十几年来其发电成本大幅下降，已具备与传统常规电源竞争的潜力。目前，国内外对于风力发电的...
复制链接

扫一扫