单边HP滤波

import pandas as pd
import numpy as np
import math
import dateutil.relativedelta
import matplotlib.pyplot as plt
import statsmodels.api as sm

#首先请自行导入一个df_original_series，格式上需要是一个仅含有一列数据的dataframe。
df_original_series.columns=['原始序列']

#%%首先定义常数矩阵At

#t:常数矩阵的维数，大于等于3的整数
#lam:HP filter的参数lambda
def At(t:int,lam):
    
    #e_t:tx1，最后一行是1，其余是0
    e_t=np.zeros(t)
    e_t[-1]=1    
    
    #I_t:t-2的单位阵
    I_t=np.identity(t-2)
    
    #Q_t:二阶差分矩阵，(t-2)xt
    Q_t=np.zeros((t-2,t)) #先设置shape
    #再通过循环设置每一行的值
    for i in range(t-2):
        Q_t[i,i],Q_t[i,i+1],Q_t[i,i+2]=1,-2,1
    
    #通过矩阵运算，计算常数阵
    # @:矩阵乘法； Matrix.T:矩阵转置； Matrix.I：矩阵求逆
    A_t=np.matrix(e_t)@(Q_t.T)@(np.linalg.inv(Q_t@(Q_t.T)+I_t/lam))@Q_t
    #结果是一个1xt的矩阵
    return A_t
    
#%%定义将序列分解为trend+cycle的函数

#输入的df：有1列待分解的数据的时间序列
#lam:HP filter的参数
def one_sided_HP_filter(df,lam):
    df_local=df.copy()
    data_series=np.array(df_local) #nx1的matrix
    length=len(df)
    
    list_cycle=[math.nan,math.nan] #t=1,2时是没有的，用math.nan填充
    for i in range(2,length): 
        #t=i+1
        sub_series=data_series[:i+1] #一共有i+1=t项
        sub_A_t=At(i+1,lam)
        cycle_t=(sub_A_t@sub_series)[0,0]
        list_cycle.append(cycle_t)
    df_local['cycle_1sHP']=list_cycle
    df_local['trend_1sHP']=df[df.columns[0]]-np.array(list_cycle)
    return df_local

#%%进行计算
df_HP=one_sided_HP_filter(df_original_series, 650)

#%%绘图，并与双侧HP滤波进行比较
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']=False

HP2s_cycle, HP2s_credit_trend = sm.tsa.filters.hpfilter(df_original_series['原始序列'],lamb=150)
df_HP['trend_2sHP']=HP2s_credit_trend

fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111)
df_HP[["原始序列", "trend_1sHP",'trend_2sHP']].plot(ax=ax, fontsize=24)
legend = ax.get_legend()
legend.prop.set_size(28)
plt.title('单侧HP滤波与双侧HP滤波的比较')
fig.savefig(path+'\图片输出\单侧HP滤波与双侧HP滤波的比较.jpg')