#include <stdint.h> #include <math.h> #include <stdio.h> #define N 200 // 最大样本量 #define p_max 5 // 最大滞后阶数 // 数据结构体 typedef struct { float data[N]; // 原始数据 float delta_y[N]; // Δy_t int n; // 实际数据长度 int p; // 滞后阶数 } ADF_Model; // OLS回归结果 typedef struct { float alpha; // 截距项α float delta; // δ系数 float gamma[p_max];// γ系数数组 float se_delta; // δ的标准误 float adf_stat; // ADF统计量 } ADF_Result; // 计算差分 Δy_t = y_t - y_{t-1} void compute_delta_y(ADF_Model *model) { for (int t = 1; t < model->n; t++) { model->delta_y[t-1] = model->data[t] - model->data[t-1]; } } // 执行ADF检验（模型2） void adf_test_model2(ADF_Model *model, ADF_Result *result) { int effective_len = model->n - 1 - model->p; // 有效样本量 float X[effective_len][model->p + 2]; // 设计矩阵 [1, y_{t-1}, Δy_{t-1}, ..., Δy_{t-p}] float y[effective_len]; // 因变量 Δy_t // 构建设计矩阵和因变量 for (int t = 0; t < effective_len; t++) { X[t][0] = 1.0; // 截距项 X[t][1] = model->data[t + model->p - 1]; // y_{t-1} y[t] = model->delta_y[t + model->p]; // Δy_t // 滞后项 Δy_{t-i} for (int i = 0; i < model->p; i++) { X[t][i + 2] = model->delta_y[t + model->p - i - 1]; } } // OLS估计（简化版矩阵运算） float XTX[(p_max + 2)][(p_max + 2)] = {0}; float XTy[(p_max + 2)] = {0}; // 计算 X'X 和 X'y for (int i = 0; i < effective_len; i++) { for (int j = 0; j < model->p + 2; j++) { for (int k = 0; k < model->p + 2; k++) { XTX[j][k] += X[i][j] * X[i][k]; } XTy[j] += X[i][j] * y[i]; } } // 矩阵求逆（使用高斯-约旦法，限于小矩阵） float inv[(p_max + 2)][(p_max + 2)]; // ... 此处省略矩阵求逆具体实现（需根据p_max优化） ... // 参数估计 beta = inv(X'X) * X'y float beta[(p_max + 2)] = {0}; for (int i = 0; i < model->p + 2; i++) { for (int j = 0; j < model->p + 2; j++) { beta[i] += inv[i][j] * XTy[j]; } } // 提取关键参数 result->alpha = beta[0]; result->delta = beta[1]; for (int i = 0; i < model->p; i++) { result->gamma[i] = beta[i + 2]; } // 计算标准误（简化方法） float residuals[effective_len]; float residual_sum = 0.0; for (int i = 0; i < effective_len; i++) { float y_pred = beta[0] + beta[1] * X[i][1]; for (int j = 0; j < model->p; j++) { y_pred += beta[j + 2] * X[i][j + 2]; } residuals[i] = y[i] - y_pred; residual_sum += residuals[i] * residuals[i]; } float sigma2 = residual_sum / (effective_len - model->p - 2); result->se_delta = sqrt(sigma2 * inv[1][1]); // ADF统计量 result->adf_stat = result->delta / result->se_delta; } // 获取p值（基于预存临界值） float get_adf_p_value(float adf_stat) { // 模型2常见临界值（样本量100） const float cv_1pct = -3.43; const float cv_5pct = -2.86; const float cv_10pct = -2.57; if (adf_stat < cv_1pct) return 0.01; else if (adf_stat < cv_5pct) return 0.05; else if (adf_stat < cv_10pct) return 0.10; else return 0.15; // 无法拒绝原假设 } int main() { ADF_Model model; ADF_Result result; // 示例数据（需替换为实际数据） float sample_data[] = { /* 时间序列数据 */ }; model.n = sizeof(sample_data)/sizeof(float); for (int i = 0; i < model.n; i++) { model.data[i] = sample_data[i]; } model.p = 2; // 滞后阶数（需根据AIC/BIC选择） compute_delta_y(&model); adf_test_model2(&model, &result); // 输出结果 printf("ADF Statistic: %.3f\n", result.adf_stat); printf("P-Value Approx: %.3f\n", get_adf_p_value(result.adf_stat)); return 0; }