HP滤波(Hodrick Prescott Filter)-java 实现

最新推荐文章于 2024-03-08 12:02:40 发布

爪洼爪洼抓光头发

最新推荐文章于 2024-03-08 12:02:40 发布

阅读量1k

点赞数

文章标签： java 开发语言

本文链接：https://blog.csdn.net/weixin_44083513/article/details/125913400

版权

HP滤波原理：

HP滤波是由Hodrick和Prescott于1980年在首先提出。这种方法被广泛地应用于对宏观经济趋势的分析研究中。HP滤波是一种时间序列在状态空间中的分析方法，相当于对波动方差的极小化。HP滤波可以看做是一个近似的高通滤波器（High-Pass Filter）。

数学公式
在这里插入图片描述

java实现方式：

依赖：

<!--矩阵处理-->
<dependency>
   <groupId>org.ejml</groupId>
   <artifactId>ejml-all</artifactId>
   <version>0.41</version>
</dependency>

依赖包api参考：
添加链接描述
 添加链接描述

 Integer lambdaPara = Integer.valueOf(formula.split(",")[0].substring(4));
 List<DeriDataPO> parse = Lists.newArrayList();
 dataMap.values().forEach(t->parse.addAll(JSONArray.parseArray(t.toString(), DeriDataPO.class)));
 List<DeriDataPO> newParse = parse.stream()
 							.sorted(Comparator.comparing(DeriDataPO::getDataDate))
 							.collect(Collectors.toList());
 List<Double> doubleList = new ArrayList<>();
 newParse.forEach(p->{doubleList.add(Double.valueOf(p.getDataValue()));});
 double[] d = new double[doubleList.size()];
 for (int i = 0; i < d.length; i++) {
     d[i] = doubleList.get(i);
 }
 Integer len = d.length;
 //创建一个(len-1) x len 的 空矩阵matrix（全部为0）
 //    SimpleMatrix(int numRows, int numCols)
 //                      行数         列数
 SimpleMatrix matrix0 = new SimpleMatrix(len - 1, len);
 //赋值
 for (int i = 0; i < len - 1; i++) {
     for (int j = 0; j < len - 1; j++) {
         if (i == j) {
            matrix0.set(i, j, -1);
         }
    }
 }
 for (int j = 1; j < len; j++) {
     for (int i = 0; i < len - 1; i++) {
         if (j - 1 == i) {
            matrix0.set(i, j, 1);
         }
    }
 }
 SimpleMatrix matrix1 = new SimpleMatrix(len - 2, len - 1);
 //赋值
 for (int i = 0; i < len - 2; i++) {
    for (int j = 0; j < len - 2; j++) {
        if (i == j) {
          matrix1.set(i, j, -1);
        }
   }
}
 for (int j = 1; j < len - 1; j++) {
     for (int i = 0; i < len - 2; i++) {
         if (j - 1 == i) {
            matrix1.set(i, j, 1);
         }
     }
 }
 //矩阵相乘
 SimpleMatrix matrix = matrix1.mult(matrix0);
 // matrix0.transpose(): A* x  A
 matrix = matrix.transpose().mult(matrix);
 //矩阵所有数乘 lambda
 //matrix.scale(lambdaPara);
 for (int i = 0; i < len; i++) {
     for (int j = 0; j < len; j++) {
         matrix.set(i , j , matrix.get(i , j)*lambdaPara);
         if (i == j){
            matrix.set(i , j , matrix.get(i , j) + 1);
         }
     }
 }
 //矩阵求逆
 SimpleMatrix matrix_invert;
 try {
     matrix_invert = matrix.invert();
 }catch (Exception e){
     throw new CommonException("计算错误");
 }
 //源数据
 SimpleMatrix sample_data_array = new SimpleMatrix(len, 1);
 for (int i = 0; i < len; i++) {
     sample_data_array.set(i, 0, d[i]);
 }
 SimpleMatrix result_trend = matrix_invert.mult(sample_data_array);
 List<Double> result = new ArrayList<>();
 for (int i = 0; i < len; i++) {
     result.add(result_trend.get(i, 0));
 }
 List<DeriDataPO> list = new ArrayList<>();
 List<DeriDataPO> dataUnify = dataUnify(JSONArray.parseArray(
 								dataMap.get(mulitfactors.get(0)).toString()));
 List<DeriDataPO> sortedByDate = dataUnify.stream()
 									.sorted(Comparator.comparing(DeriDataPO::getDataDate))
 									.collect(Collectors.toList());
for (int i = 0; i < sortedByDate.size(); i++) {
    DeriDataPO deriDataPO = new DeriDataPO();
    deriDataPO.setIndexCode(formula);
    deriDataPO.setDataDate(sortedByDate.get(i).getDataDate());
    deriDataPO.setDataValue(result.get(i).toString());
    list.add(deriDataPO);
}
return list;

python实现方式：

###输入的时间序列
sample_data=[4137.8,4109.4,4197.4,4153.2,4082.4,3952.6,4017.0,
4028.4,4015.0,3983.0,4067.2,4178.0,4277.6,4606.2,
4879.4,4855.2, 4698.6,4422.4,4303.25,3995.2,3883.8,
4003.6,3983.2,4016.6,4124.0,4472.0,4483.6,4137.8,4109.4,4197.4,4153.2,4082.4,3952.6,4017.0,
4028.4,4015.0,3983.0,4067.2,4178.0,4277.6,4606.2,
4879.4,4855.2, 4698.6,4422.4,4303.25,3995.2,3883.8,
4003.6,3983.2,4016.6,4124.0,4472.0,4483.6,
4137.8,4109.4,4197.4,4153.2,4082.4,3952.6,4017.0,
4028.4,4015.0,3983.0,4067.2,4178.0,4277.6,4606.2,
4879.4,4855.2, 4698.6,4422.4,4303.25,3995.2,3883.8,
4003.6,3983.2,4016.6,4124.0,4472.0,4483.6]
lambda_para=5          ##超参数，前端输入，平滑系数，越大越平滑，这个暴露给外面作为参数
N=len(sample_data) #输入数据的长度，float类型且不为空
F_para1=np.zeros([N-1,N])  ##构建一阶导数系数
for i in range(0,N-1):
for j in range(0,N-1):
if i==j:
F_para1[i,j]=-1
for j in range(1,N):
for i in range(0,N-1):
if j-1==i:
F_para1[i,j]=1
F_para2=np.zeros([N-2,N-1])
for i in range(0,N-2):
for j in range(0,N-2):
if i==j:
F_para2[i,j]=-1
for j in range(1,N-1):
for i in range(0,N-2):
if j-1==i:
F_para2[i,j]=1
F_para=np.matmul(F_para2,F_para1) ##矩阵相乘
F_para=np.matmul(F_para.T,F_para) ##矩阵相乘
for i in range(N):
for j in range(N):
F_para[i,j]=F_para[i,j]* lambda_para
if i==j:
F_para[i,j]=F_para[i,j]+1
F_para_inverse=np.linalg.inv(F_para) ##矩阵求逆
sample_data_array=np.zeros([N,1])
for i in range(N):
sample_data_array[i,0]=sample_data[i]
result_trend=np.matmul(F_para_inverse,sample_data_array)   #矩阵相乘        #趋势项
result_trend=result_trend[:,0]                                     #趋势项目
result_resid=sample_data_array[:,0]-result_trend                   #残差项
plt.figure()
plt.plot(sample_data_array[:,0], label='sorce_data')
plt.plot(result_trend, label='trend')
plt.plot(result_resid, label='resid')
plt.legend()