疫情期間個類別股票的表現以及因子模型的失效

Posted on by Jacky Hsueh

2020開始的COVID-19對於金融市場其實是一個巨大的衝擊，而疫情期間其實許多很多現象都不太一樣，這一篇就記錄一下2020-2022的資料以及先前的資料預測不太一樣的地方。

基本上資料展項兩個有趣的地方:1.產業之間的共變異數改變，2.價值投資的失靈。我個人覺得挺有趣的，就寫一下當作紀錄。

1.產業之間共變異數改變

在疫情發生以前，如果我們用各種ETF的價格來當作個產業的報酬，把產業分為:通訊(VOX), 非必要性消費(VCR), 必要性消費(VDC), 能源(VDE), 金融(VFH), 健康(VHT), 工業(VIS), 科技(VGT), 原物料(VAW), 房地產(VNQ), 水電(VPU)。並且比較個產業相較於SP500指數(SPY)的關係，我們可以得到下面的圖一:

圖一:各產業對於彼此之間報酬的相關係數，由下方從左到右分別為:
能源、健康、科技、金融、原物料、必要性消費、工業、SP500、通訊、非必要性消費、房地產、水電

雖然說SP500本來就是各產業的組成，但有些產業其實跟大盤的連動比較大(例如科技、工業)，有些產業則是沒那麼明顯(例如房地產、水電)。而下面圖二則是各種產業以及SP500在2010-2020期間每一季的利潤相關係數的變化，而我們可以看到在有些年分，其實原本高相關係數的反而會跟SP500的相關係數沒有那麼多的關係

圖二:各產業與SP500的相關係數

圖三可以比較房地產以及科技這兩種股票歷年來跟SP500報酬的關係。事實上，我們可以看到SP500的報酬在某些年分與房地產還有科技類股票的關聯性其實不大，例如在2018年左右，房地產以及科技類股票與SP500的報酬其實只有0.2以及0.6的。

圖三:房地產(紫色)以及科技業(藍色)相對於SP500的關係

那麼疫情期間有什麼改變呢?下圖的圖四比較了疫情以前以及疫情以後的改變，上面的圖是疫情發生以前各產業對於SP500的相關係數，我們可以看到其實原本最高的是工業，其次是必要性消費，接著才是科技。但是在疫情發生以後，順序發生了改變，工業以及SP500的關係反而沒那麼重要，倒是通訊產業以及SP500的相關性從原本只有0.7左右來到了將近0.9。而總體而言，原本只有6個產業與SP500的報酬是高度相關的(>0.8)，在疫情期間反而變成9種產業。

原因或許在於疫情期間中央銀行大量放水，導致資金過多，使得股是裡面所有產業都幾乎是能賺錢的，導致了各產業都隨著SP500一起變動。而值得注意的是在疫情期間反而能源股的報酬以及SP500的相關係數不大。

圖四:疫情發生前各產業與SP500的相關係數(上)以及疫情發生後各產業與SP500的相關係數(下)

2.失靈的因子模型

接著，我跑了因子模型分別針對各產業去作兩階段回歸。所謂的因子模型，出自於Fama以及French兩位金融學家的創舉，基本上可以看做是量化版本的價值投資，他們找出五種可以去分析公司股票的方式，然後歸法為五種因子，而根據這五種因子去解釋為什麼有些公司或者產業可以產出超額報酬。這五種因子分別為:Mkt-RF(市場報酬率，定義為市場上所有股票的報酬扣掉國庫債報酬)、SMB(大小報酬利差，定義為投資組合裡面九個最小的報酬資產扣掉九個最大報酬資產的利差)、 HML(高低利差，定義為兩種價值股的平均報酬以及兩種成長股報酬的利差)、 RMW(穩利與不穩定利差，定義為兩種穩定報酬組合的平均以及兩種高風險報酬組合平均的利差)、CMA (風險與激進利差，定義為兩種保守的投資組合的平均以及兩種激進投資組合平均的利差)。

基本上，量化基金例如AQR就是使用這種方式去尋找值得投資的價值股票的代表，而我們已經看到了各產業在不同期間其實相對於SP500會有不同的變異數，那麼我想問的就是如何透過這種方式來分析哪些產業相較於SP500能產生超額報酬呢?換句話說，我想確定這種投資方法在疫情發生以前以及疫情發生以後的效果如何。

基本上我使用的方法被稱之為Fama-MacBeth，基本的概念就是我先回歸每一種資產的抽額報酬以及他的因子，例如我們總共有17個資產，而橫跨了T個時間點，那針對每一個資產我都可以產生T個回歸，而因為我有5個因子，加上常數項的話則是有6個因子要跑，所以我可以跑一個T*6的矩陣，而每個資產可以得到6個coefficient。

然而，我們需要知道各種資產的以及因子的風險溢酬，所以我還要再跑一次回歸，這一次我直接把所有資產放進來跑回歸，因此我要跑的變數會是17*6的矩陣代表所有資產的差額報酬，而獨立變數則會是剛剛得到的17*6矩陣的coefficient，而得到的新的coefficient才會是能夠有效解釋各種因子風險溢酬的回歸。

之所以要跑這樣兩段式回歸的原因在於如果你單跑第二段的回歸，你會把不同資產的變異數以及誤差重新再算一次，所以我們要先分開跑，接著再一起跑，如此一來就可以消除多項資產之間變異數的問題，這方法就是著名的Fama-MacBeth。

而下圖五以及表一則是疫情發生以前各種因子他的比重以及回歸表:

表一:2010-2020因子投資回歸
圖五:2010-2020因子投資比重

我們首先來介紹一下怎麼讀統計表，可以直接去看那個P-Value，也就是拒絕這項因子可以解釋超額報酬的機率，可以看到市場報酬(Mkt-RF)、高低利差(HML)以及保守與激進利差(CMA)其實是可以解釋的，因為他們拒絕的機率其實只有5%左右，而這三種的比重也的確占很重。換句話說，在2010-2020這段期間裡面，我們其實可以套過這樣的因子投資來事先找出可以獲利的公司以及產業其實理論上是沒問題的，甚至可以事先知道自己的投資策略是否夠穩。

那麼2020疫情爆發以後呢?圖六以及表二呈現了回歸的結果:

表二:疫情期間因子投資回歸
圖六:疫情期間因子投資比重

我們可以看到那個P-Value基本上幾乎都是慘不忍睹，基本上都是統計上面沒啥意義的大，換句話說，基本上因子投資在這段期間幾乎無法用來有效的預測哪些產業可以產生超額報酬。

即便撇開統計顯著不談，我們來比較圖五以及圖六，會發現有些coefficient原本是負的(CMA，RMW)在疫情期間變成正的，而原本比重不大的SMB突然比重拉大，因此相對而言HML因子的比重是相對被拉小的。也就是說，如果你用2010-2020很穩的價值投資策略去操盤2020-2022，那會死很慘….因為因子的方向完全不對。

結論:過去策略應該不斷回測

前面提到大量使用因子投資模型的量化公司AQR在2020年發了一篇”Is (Systematic) Value Investing Dead?“的文章講樹價值投資是不是死了，而至少我們從2020-2022的資料來看，很難說這些價值因子可以跟2020以前一樣獲得穩定的報酬。

我想，剛好自己研究這個的啟示大概就是:”過去的策略需要不斷地回測”，當總體經濟產生改變的時候，我們也要改變自己的想法。有趣的研究方向或使是去問說，那是什麼造成這種風向大轉彎的呢?是因為不同產業與SP500的變異數改變了?還是因為低利率?或者是因為其他高風險產品的出現改變了人們對於利潤的定錨?

我不知道，但我會很想知道，而我也在朝這麼方向研究當中，這一篇就當作紀錄一下，未來或許有機會可以回來看。

附錄:程式碼

from pprint import pprint
import numpy as np
import pandas as pd
import pyfolio as pf
import quantstats as qs
from datetime import date
from pandas_datareader.famafrench import get_available_datasets
from scipy.stats import spearmanr, pearsonr
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from statsmodels.api import OLS, add_constant
from pathlib import Path
from linearmodels.asset_pricing import TradedFactorModel, LinearFactorModel, LinearFactorModelGMM
import seaborn as sns 
from pandas_datareader import data as web
import warnings
warnings.filterwarnings('ignore')
pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', 500)

ff_portfolio = '17_Industry_Portfolios'
ff_portfolio_data = web.DataReader(ff_portfolio, 'famafrench', start=start, end=end)[0]
ff_portfolio_data = ff_portfolio_data.sub(ff_factor_data.RF, axis=0)
ff_portfolio_data.info()

symbols = ['VOX', 'VCR', 'VDC', 'VDE', 'VFH', 'VHT', 'VIS', 'VGT', 'VAW', 'VNQ', 'VPU']
secs = ['COMM', 'CONSUMER DISC', 'CONSUMER ST', 'ENERGY', 'FINANCIALS', 'HEALTH', 'INDUSTRIALS', 
           'TECHNOLOGY', 'MATIREALS', 'REAL ESTATE', 'UTILITIES']

def get_symbols(symbols,data_source,ohlc,begin_date=None,end_date=None):
    out = []
    new_symbols = []
    
    for symbol in symbols:
        df = web.DataReader(symbol, data_source,begin_date, end_date)\
        [['High','Low','Open','Close','Volume','Adj Close']]
        new_symbols.append(symbol) 
        out.append(df[ohlc].astype('float'))
        data = pd.concat(out, axis = 1)
        data.columns = new_symbols
        
    return data

prices = get_symbols(symbols,data_source='yahoo',ohlc='Close',\
                     begin_date=start,end_date=end)

SPY = web.DataReader('SPY', 'yahoo', start, end)\
      [['High','Low','Open','Close','Volume','Adj Close']]

prices.columns = secs
prices['SPY'] = SPY.Close.values
returns0 = prices.pct_change()
returns0 = returns0.dropna(how='all').dropna(axis=1)

上面為準備資料的，下面產出圖一到四以及其他圖片

pricesx=returns0
ricesx = returns0.copy()

n_secs = len(returns0.columns)
colors = cm.rainbow(np.linspace(0, 1, n_secs))
returns0.cumsum().plot(color=colors, figsize=(12, 6))# Normalize Prices 
plt.title('Cummulative Returns Series')
plt.xlabel('Date')
plt.ylabel('Returns (%)')
plt.grid(b=None, which=u'major', axis=u'both')
plt.legend(bbox_to_anchor=(1.01, 1.1), loc='upper left', ncol=1)
plt.show();


pricesx.hist(bins=20, figsize=(10,10))
plt.tight_layout()
plt.show()

plt.figure(figsize=(12,6))
g1 = sns.boxplot(data=pricesx)
g1.set_xticklabels(g1.get_xticklabels(),rotation=90)
plt.title('Returns')
plt.tight_layout()
plt.show()

g = sns.clustermap(pricesx.corr(), annot=True)
ax = g.ax_heatmap
bottom, top = ax.get_ylim()
ax.set_ylim(bottom + 0.5, top - 0.5)
plt.title('Pairwise Correlations')
plt.show()

pricesx.drop('SPY', axis=1).rolling(63).corr(pricesx.SPY).dropna().plot(color=colors, figsize=(12, 6))
plt.legend(bbox_to_anchor=(1.01, 1.1), loc='upper left', ncol=1)
plt.title('Rolling Quarterly Correlation to SPY')
plt.show()

pricesx.drop('SPY',axis=1).corrwith(pricesx.SPY).sort_values(ascending=True).plot(kind='barh', color=colors, figsize=(12, 6))
plt.title('Correlation to SPY')
plt.show()

pricesx.rolling(63).var().dropna().plot(color=colors, figsize=(12, 6))
plt.legend(bbox_to_anchor=(1.01, 1.1), loc='upper left', ncol=1)
plt.title('Rolling Quarterly Variance')
plt.show()

plt.figure(figsize=(12,6))
g2 = sns.boxplot(data=pricesx.rolling(63).std().dropna())
g2.set_xticklabels(g2.get_xticklabels(),rotation=90)
plt.title('Rolling Quarterly Returns Volatility')
plt.tight_layout()
plt.show()


pricesx.drop('SPY', axis=1).rolling(63).corr(pricesx.SPY).dropna().plot(color=colors, figsize=(12, 6))
plt.legend(bbox_to_anchor=(1.01, 1.1), loc='upper left', ncol=1)
plt.title('Rolling Quarterly Correlation to SPY')
plt.show()

pricesx.drop('SPY',axis=1).corrwith(pricesx.SPY).sort_values(ascending=True).plot(kind='barh', color=colors, figsize=(12, 6))
plt.title('Correlation to SPY')
plt.show()

pricesx.rolling(63).var().dropna().plot(color=colors, figsize=(12, 6))
plt.legend(bbox_to_anchor=(1.01, 1.1), loc='upper left', ncol=1)
plt.title('Rolling Quarterly Variance')
plt.show()

plt.figure(figsize=(12,6))
g2 = sns.boxplot(data=pricesx.rolling(63).std().dropna())
g2.set_xticklabels(g2.get_xticklabels(),rotation=90)
plt.title('Rolling Quarterly Returns Volatility')
plt.tight_layout()
plt.show()

pricesx.get(['REAL ESTATE','TECHNOLOGY']).rolling(63).corr(pricesx.SPY).dropna().plot(color=colors, figsize=(12, 6))
plt.title('Rolling Quarterly Correlation to SPY')
plt.show()

圖1-4結束，下面是產出Fama-MacBeth的

def r2(x, y):
    return pearsonr(x, y)[0] ** 2

tmp_fts = pricesx.copy()
returns = prices.resample('M').last().pct_change().mul(100).to_period('M')
returns = returns.dropna(how='all').dropna(axis=1)
returns.info()

ff_factor_data = ff_factor_data.loc[returns.index]
ff_portfolio_data = ff_portfolio_data.loc[returns.index]
ff_factor_data.describe()

excess_returns = returns.sub(ff_factor_data.RF, axis=0)
excess_returns.info()
excess_returns = excess_returns.clip(lower=np.percentile(excess_returns, 1),
                                     upper=np.percentile(excess_returns, 99))
ff_portfolio_data.info()

from statsmodels.api import OLS, add_constant

betas = []
for industry in ff_portfolio_data:
    step1 = OLS(endog=ff_portfolio_data.loc[ff_factor_data.index, industry], 
                exog=add_constant(ff_factor_data)).fit()
    betas.append(step1.params.drop('const'))

betas = pd.DataFrame(betas, 
                     columns=ff_factor_data.columns, 
                     index=ff_portfolio_data.columns)
betas.info()

lambdas = []
for period in ff_portfolio_data.index:
    step2 = OLS(endog=ff_portfolio_data.loc[period, betas.index], 
                exog=betas).fit()
    lambdas.append(step2.params)

lambdas = pd.DataFrame(lambdas, 
                       index=ff_portfolio_data.index,
                       columns=betas.columns.tolist())
lambdas.info()

ax1 = plt.subplot2grid((1, 3), (0, 0))
ax2 = plt.subplot2grid((1, 3), (0, 1), colspan=2)
ax2.margins(0.01)
lambdas.mean().plot.barh(ax=ax1)
lambdas0 = lambdas.rolling(6).mean().dropna()
lambdas0.plot(lw=2, figsize=(17,8), ax=ax2)
ax2.legend(bbox_to_anchor=(1.025, 1.05))
plt.show()

lambdas.rolling(12).mean().dropna().plot(lw=2, figsize=(14,20), subplots=True, sharey=True, sharex=True)
plt.show()

上面是產生圖5,6，下面是產生表1,2:

from linearmodels.asset_pricing import TradedFactorModel, LinearFactorModel, LinearFactorModelGMM

mod = LinearFactorModel(portfolios=ff_portfolio_data, 
                        factors=ff_factor_data)
res = mod.fit(cov_type='robust')
print(res.summary)

發佈留言

發佈留言必須填寫的電子郵件地址不會公開。 必填欄位標示為 *