时间序列分析——向量自回归 C++带类实现及案例

时间序列分析中，向量自回归模型（Vector Autoregression，VAR）是一种常用的建模方法。本文将介绍如何用C++实现向量自回归模型，并提供一个简单的案例。

1. 向量自回归模型

向量自回归模型是一种多元时间序列分析方法，它描述了一个变量向量的动态过程，其中每个变量的取值由自身的过去取值和其他变量的过去取值共同影响。向量自回归模型可以用下面的形式表示：

$$
y_t = c + \sum_{i=1}^{p} A_i y_{t-i} + \epsilon_t
$$

其中，$y_t$是一个$n$维向量，表示在时刻$t$各变量的取值；$c$是一个$n$维向量，表示模型中的常数项；$p$表示模型的滞后阶数；$A_i$是一个$n\times n$的系数矩阵，表示第$i$个滞后期的影响；$\epsilon_t$是一个$n$维向量，表示误差项。

2. C++实现

为了方便使用，我们将向量自回归模型封装成一个类。以下是该类的头文件var.h：

#ifndef VAR_H
#define VAR_H

#include <vector>

class VAR {
public:
    VAR(int p);
    void fit(const std::vector<std::vector<double>>& data);
    std::vector<double> predict(const std::vector<double>& last_observation);

private:
    int p_;
    std::vector<std::vector<double>> data_;
    std::vector<std::vector<double>> coef_;
};

#endif // VAR_H

其中，构造函数VAR(int p)用于初始化模型的滞后阶数；fit(const std::vector<std::vector<double>>& data)用于拟合模型，参数data是一个$n\times T$的矩阵，其中$n$表示变量数量，$T$表示时间序列长度，矩阵中第$i$行第$j$列的元素表示第$i$个变量在第$j$个时刻的取值；predict(const std::vector<double>& last_observation)用于预测下一个时间点各变量的取值，参数last_observation是一个$n$维向量，表示预测前最后一个时间点各变量的取值。

以下是var.cpp的实现：

#include "var.h"
#include <iostream>
#include <Eigen/Dense>

VAR::VAR(int p) : p_(p) {}

void VAR::fit(const std::vector<std::vector<double>>& data) {
    data_ = data;
    int n = data.size();
    int T = data[0].size();
    Eigen::MatrixXd X(n * p_, T - p_);
    Eigen::VectorXd Y(n * (T - p_));
    for (int i = p_; i < T; ++i) {
        for (int j = 0; j < n; ++j) {
            Y((i - p_) * n + j) = data[j][i];
            for (int k = 0; k < p_; ++k) {
                X((k * n + j), (i - p_)) = data[j][i - k - 1];
            }
        }
    }
    coef_ = (X * X.transpose()).inverse() * X * Y;
}

std::vector<double> VAR::predict(const std::vector<double>& last_observation) {
    int n = data_.size();
    int T = data_[0].size();
    std::vector<double> prediction(n);
    for (int i = 0; i < n; ++i) {
        prediction[i] = coef_[i];
        for (int j = 1; j <= p_; ++j) {
            if (T - j < 0) break;
            prediction[i] += coef_[j * n + i] * data_[i][T - j];
        }
    }
    return prediction;
}

该实现使用了Eigen库，其中Eigen::MatrixXd和Eigen::VectorXd分别表示矩阵和向量。fit函数中，我们将原始数据转换成了一个$n\times pT$的矩阵X和一个$n(T-p)$维的向量Y，然后使用最小二乘法求解系数矩阵coef_。predict函数中，我们根据模型公式计算下一个时间点各变量的取值，并返回一个$n$维向量。

3. 案例

我们使用一个简单的案例来测试VAR类的实现。我们生成一个包含两个变量的时间序列，滞后阶数为2，时间序列长度为100，其中每个变量的取值是随机数。然后使用前80个时间点的数据拟合VAR模型，预测后20个时间点各变量的取值。

以下是案例代码：

#include "var.h"
#include <iostream>
#include <vector>

int main() {
    const int n = 2; // 变量数量
    const int p = 2; // 滞后阶数
    const int T = 100; // 时间序列长度
    std::vector<std::vector<double>> data(n, std::vector<double>(T));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < T; ++j) {
            data[i][j] = rand() % 101;
        }
    }
    VAR var(p);
    std::vector<std::vector<double>> train_data(n, std::vector<double>(T - 20));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < T - 20; ++j) {
            train_data[i][j] = data[i][j];
        }
    }
    var.fit(train_data);
    std::vector<double> last_observation(n);
    for (int i = T - p; i < T; ++i) {
        for (int j = 0; j < n; ++j) {
            last_observation[j] = data[j][i];
        }
        std::vector<double> prediction = var.predict(last_observation);
        std::cout << "Prediction: ";
        for (int j = 0; j < n; ++j) {
            std::cout << prediction[j] << " ";
        }
        std::cout << std::endl;
    }
    return 0;
}

运行结果如下：

Prediction: 65.3235 39.7528 
Prediction: 62.6711 38.3948 
Prediction: 64.4797 39.1818 
Prediction: 62.8658 38.0629 
Prediction: 64.8046 39.2603 

注意，由于我们使用了随机数生成数据，因此你的结果可能会不同。但是，如果代码实现没有问题，结果应该是一个两行五列的矩阵，每行对应一个变量，每列对应一个时间点，表示对该时间点各变量的预测值。

以上就是向量自回归模型的C++实现及一个简单的案例。