分位数VAR模型¶

import pandas as pd
import statsmodels.regression.quantile_regression as qr
import statsmodels.api as sm

data = pd.read_excel('../数据/上证指数与沪深300.xlsx')
Y = data['hs300']
X = data['sz']

def lag_list(Y, X, p=1, q=1, Yname='', Xname='', exogenous=[]):
    '''
    待估计方程：y = c + y(-1) +....+y(-p) + x(-1) + ... + x(-q) + exogenous
    获取自回归分布滞后模型的估计向量

    Parameters
    ----------
    Y : 被估计变量
    X : 估计变量
    p : ADL 模型 Y 的滞后阶数，标量默认为1
    q : ADL 模型 X 的滞后阶数，标量默认为1
    Yname : 被响应变量名
    Xname : 响应变量名
    exogenous : 外生变量

    Returns
    -------
    ADLy : ADL 模型被解释变量
    ADLx : ADL 模型解释变量

    '''
    if not Yname:
        Yname = 'y'
    if not Xname:
        Xname = 'x'
        
    ADLx = pd.DataFrame()
    T = len(Y)
    ADLy = pd.Series(list(Y[max(p, q):T]), name=Yname)
    for i in range(1, p+1):
        name = f'{Yname}_{i}'
        ADLx[name] = list(Y[max(p, q)-i:T-i])
    for i in range(1, q+1):
        name = f'{Xname}_{i}'
        ADLx[name] = list(X[max(p, q)-i:T-i])
    
    # 增加控制变量
    if type(exogenous) == pd.Series:
        ADLx[exogenous.name] = list(exogenous[:0-max(p, q)])
    elif type(exogenous) == pd.DataFrame:
        for name in exogenous.columns:
            ADLx[name] = list(exogenous[name][:0-max(p, q)])
    
    # 增加常数项
    ADLx = sm.add_constant(ADLx)
    return ADLy, ADLx


def qvar(Y, X, P, Q, Yname='', Xname='', exogenous=[]):
    '''
    待估计方程：y = c + y(-1) +....+y(-p) + x(-1) + ... + x(-q) + exogenous
               x = c + y(-1) +....+y(-p) + x(-1) + ... + x(-q) + exogenous
    估计 QVAR 模型

    Parameters
    ----------
    Y : 被估计变量
    X : 估计变量
    P : QVAR 模型的滞后阶数
    Q : 分位点
    Yname : 被响应变量名
    Xname : 响应变量名
    exogenous : 外生变量
    
    Returns
    -------
    res1, res2 : 两模型估计结果
    
    '''
    ADLy, ADLx = lag_list(Y, X ,P, P, Yname, Xname, exogenous)
    mod = qr.QuantReg(ADLy, ADLx)
    res1 = mod.fit(Q)
    
    ADLy, ADLx = lag_list(X, Y ,P, P, Xname, Yname, exogenous)
    mod = qr.QuantReg(ADLy, ADLx)
    res2 = mod.fit(Q)
    
    return res1, res2

res1, res2 = qvar(Y, X, 2, 0.3, 'hs300', 'sz')
beta1 = res1.params
beta2 = res2.params
pd.DataFrame([beta1, beta2], ['hs300', 'sz'])

import numpy as np

def OIRF(res1, res2, p):
    '''
    估计脉冲响应函数

    Parameters
    ----------
    res1, res2 : 两模型估计结果
    p : 滞后阶数
    
    Returns
    -------
    OIRF
    
    '''
    pass
p = 2
resid = pd.DataFrame([res1.resid, res2.resid]).T   # 残差序列
# 正交化分解，估计P2矩阵
a = resid.T @ resid
a=a/(len()-2*p-1);
P2 = chol(a, 'lower');