import numpy as np from arch import arch_model def run_univariate_garch(ret): model_params = {} T,N = ret.shape D= np.zeros((T,N)) for r_i in ret: am = arch_model (ret[r_i]) model_params[r_i] = am.fit(disp='off') D[:,ret.columns.get_loc(r_i)] =model_params[r_i].conditional_volatility return model_params, D # calculate R_t and Q as in quation (1) return Rt and vech(R_t) - lower part of matrix as stacked vector def dcceq(phi,ret,D_inv,init, Eps = None): a, b = phi T,N = ret.shape if Eps is None: Eps = np.zeros((T,N)) for t in range(T): Eps[t,:] = np.matmul(np.diagflat(D_inv[t,:]), ret[t]) S= 1/T*np.matmul(Eps.T,Eps) if min(a,b)<0 or max(a,b)>1 or a+b > .999999: a = .9999 - b Qt = np.zeros((N, N ,T)) Qt[:,:,0] = np.cov(Eps[:init,:].T) Eps_shifted = np.insert(Eps, 0, Eps[0,:], axis=0) Eps_shifted = np.delete(Eps_shifted,obj = -1, axis =0) Rt = np.zeros((N, N ,T)) Rt[:,:,0] = np.corrcoef(Eps[:init,:].T) for n in range(0,N): for m in range(0,N): Qt[n,m,:] = S[n,m] * (1-a-b) Qt[n,m,:] = Qt[n,m,:] + a * Eps_shifted[:,n]*Eps_shifted[:,m] Qt[:,:,0] = np.cov(Eps[:init,:].T) for t in range(1,T): Qt[:,:,t] = Qt[:,:,t] + b* Qt[:,:,t-1] for n in range(0,N): for m in range(0,N): Rt[n,m,:] = np.divide(Qt[n,m,:],np.sqrt(Qt[n,n,:]*Qt[m,m,:])) return Rt def loglike_dcc_svd(phi,ret,D_inv,init, SVD= False, verbose= False): T,N = ret.shape llf = np.zeros((T,1)) Eps = np.zeros((T,N)) for t in range(T): Eps[t,:] = np.matmul(np.diagflat(D_inv[t,:]), ret[t]) Rt = dcceq(phi,ret,D_inv,init,Eps) U = np.zeros((N, N ,T)) s = np.zeros((N ,T)) Vh = np.zeros((N, N ,T)) if SVD == True: for j in range(0,T): U[:,:,j],s[:,j],Vh[:,:,j]= np.linalg.svd(Rt[:,:,j], full_matrices=True, compute_uv=True, hermitian=True) if SVD == True: for i in range(0,T): rV = np.matmul(np.array(Eps[i,:], ndmin=2) , Vh[:,:,i].T) Ur = np.matmul(U[:,:,i].T, np.array(Eps[i,:], ndmin=2).T) rVS = np.divide(rV, s[:,i]) det = np.log(s[:,i]).sum() mult = np.matmul(rVS, Ur)- np.matmul(np.array(Eps[i,:], ndmin=2), np.array(Eps[i,:], ndmin=2).T) #Make sure this makes sense llf[i] =det + mult llf = np.sum(llf) else: for i in range(0,T): det = np.log(np.linalg.det(Rt[:,:,i])) mult = np.matmul(np.matmul(Eps[i,:].T , (np.linalg.inv(Rt[:,:,i]) - np.eye(N))) ,Eps[i,:]) llf[i] =(det+ mult) llf = np.sum(llf) if verbose: print(phi,llf) return llf def composite_pair_loglike_dcc(phi,ret,D_inv,init): T,N = np.shape(ret) Rj = np.zeros((N, N ,T)) llfj = np.zeros((T,1)) Eps = np.zeros((T,N)) for t in range(T): Eps[t,:] = D_inv[t,:]* ret[t] Eps_sq = Eps**2 Rj = dcceq(phi,ret,D_inv,init, Eps) Rj_inv = np.zeros((N,N,T)) det = np.zeros((T)) det[:] = Rj[0,0,:]* Rj[1,1,:]- Rj[0,1,:]* Rj[1,0,:] Rj_inv[0,0,:] = Rj[1,1,:] /det[:] Rj_inv[1,1,:] = Rj[0,0,:] /det[:] Rj_inv[0,1,:] = -Rj[0,1,:] /det[:] Rj_inv[1,0,:] = -Rj[1,0,:] /det[:] Eps_sq_col_sum = Eps_sq.sum(axis = 1) mult = np.einsum('ij, jki, ik ->', Eps, Rj_inv, Eps) mult_total = mult - np.sum(Eps_sq_col_sum) llfj = np.sum(np.log(det))+np.sum(mult_total) return llfj def composite_loglike(phi,ret,D_inv,init, endpair= None, verbose = False): T,N = np.shape(ret) llf = np.zeros((int(N*(N-1)/2),1)) i=0 if endpair is not None: endpair = min(N,endpair) endpair = max(2,endpair) else: endpair = N for m in range(1,endpair): for n in range(0,N): if m+n < N: retj =ret[:,[n,m+n]] D_invj = D_inv[:,[n,m+n]] llf[i]= composite_pair_loglike_dcc(phi,retj,D_invj,init) i+=1 llf= np.sum(llf) if verbose: print(phi,llf) if llf == float("-inf"): #trick to keep alhorithom going even when singular matrix is obtained in some interation llf=0 return llf ---------------------------------------------------------------------------------------------------- import numpy as np from arch import arch_model def run_univariate_garch(ret): model_params = {} T,N = ret.shape D= np.zeros((T,N)) for r_i in ret: am = arch_model (ret[r_i]) model_params[r_i] = am.fit(disp='off') D[:,ret.columns.get_loc(r_i)] =model_params[r_i].conditional_volatility return model_params, D def ccceq(ret, D_inv,init): T,N = ret.shape Eps = np.zeros((T,N)) for t in range(T): Eps[t,:] = np.matmul(np.diagflat(D_inv[t,:]), ret[t]) Rt = np.zeros((N, N ,T)) for j in range(0,min(init,N)): Rt[:,:,j] = np.cov(Eps[:init,:].T) for j in range(min(init,N),T): Rt[:,:,j] = np.cov(Eps[:(j-1),:].T) return Rt ----------------------------------------------------------------------------------------- import numpy as np def ewma_est_H(a,ret,init): T,N = np.shape(ret) Ht = np.zeros((N, N ,T)) Ht[:,:,0] = np.cov(ret[:init,:].T) for j in range(1,T): Ht[:,:,j] = a * np.matmul(np.array(ret[j-1,:], ndmin=2).T,np.array(ret[j-1,:], ndmin=2)) Ht[:,:,j] = Ht[:,:,j] + (1-a) * Ht[:,:,j-1] return Ht def eq_ewma_svd(a,ret,init, SVD=False): T,N = np.shape(ret) Ht = np.zeros((N, N ,T)) U = np.zeros((N, N ,T)) s = np.zeros((N ,T)) V = np.zeros((N, N ,T)) ret_shifted = np.insert(ret, 0, ret[0,:], axis=0) ret_shifted = np.delete(ret_shifted,obj = -1, axis =0) Ht[:,:,0] = np.cov(ret[:init,:].T) if SVD == True: U[:,:,0],s[:,0],V[:,:,0] = np.linalg.svd(Ht[:,:,0], full_matrices=True, compute_uv=True, hermitian=True) for n in range(0,N): for m in range(0,N): Ht[n,m,:] = a * ret_shifted[:,n]*ret_shifted[:,m] Ht[:,:,0] = np.cov(ret[:init,:].T) for t in range(1,T): Ht[:,:,t] = Ht[:,:,t] + (1-a)* Ht[:,:,t-1] if SVD == True: U[:,:,t],s[:,t],V[:,:,t] = np.linalg.svd(Ht[:,:,t], full_matrices=True, compute_uv=True, hermitian=True) if SVD == True: return U,s,V else: return Ht def loglike_ewma_svd(a,ret,init, SVD = False, verbose = False): T,N = np.shape(ret) llf = np.zeros((T,1)) outs =[] if SVD == True: U,s,Vh = eq_ewma_svd(a,ret,init,SVD) else: Ht= eq_ewma_svd(a,ret,init, SVD) if SVD == True: for i in range(0,T): rV = np.matmul(np.array(ret[i,:], ndmin=2) , Vh[:,:,i].T) Ur = np.matmul(U[:,:,i].T, np.array(ret[i,:], ndmin=2).T) rVS = np.divide(rV, s[:,i]) llf[i] = np.matmul(rVS, Ur) llf = np.sum(llf) + s.sum() else: for i in range(0,T): det = np.log(np.linalg.det(Ht[:,:,i])) mult = np.matmul(np.matmul(np.array(ret[i,:], ndmin=2) , (np.linalg.inv(Ht[:,:,i]))) ,np.array(ret[i,:], ndmin=2).T) llf[i] =(det+ mult) llf = np.sum(llf) if verbose: print(a,llf) if llf == float("-inf"): #trick to keep alhorithom going even when singular matrix is obtained in some interation llf=0 return llf def composite_pair_loglike_ewma(a,ret,init): T,N = np.shape(ret) Hj = np.zeros((N, N ,T)) llfj = np.zeros((T,1)) Hj = eq_ewma_svd(a,ret,init, False) Hj_inv = np.zeros((N,N,T)) det = np.zeros((T)) det[:] = Hj[0,0,:]* Hj[1,1,:]- Hj[0,1,:]* Hj[1,0,:] Hj_inv[0,0,:] = Hj[1,1,:] /det[:] Hj_inv[1,1,:] = Hj[0,0,:] /det[:] Hj_inv[0,1,:] = -Hj[0,1,:] /det[:] Hj_inv[1,0,:] = -Hj[1,0,:] /det[:] mult = np.einsum('ij, jki, ik ->', ret, Hj_inv, ret) llfj = np.sum(np.log(det))+np.sum(mult) return llfj def composite_loglike(a,ret,init, endpair= None, verbose =False): T,N = np.shape(ret) llf = np.zeros((int(N*(N-1)/2),1)) i=0 if endpair is not None: endpair = min(N,endpair) endpair = max(2,endpair) else: endpair = N for m in range(1,endpair): for n in range(0,N): if m+n < N: retj =ret[:,[n,m+n]] llf[i]= composite_pair_loglike_ewma(a,retj,init) i+=1 llf= np.sum(llf) if verbose: print(a,llf) return llf ------------------------------------------------------------------------------------------------ import numpy as np def svech_est_H(theta,ret,init,intercept_m= None): T,N = np.shape(ret) a,b = theta Ht = np.zeros((N, N ,T)) if intercept_m is None: intercept_m = np.zeros((N,N)) for t in range(T): intercept_m= intercept_m +np.matmul(np.array(ret[t,:], ndmin=2).T,np.array(ret[t,:], ndmin=2)) intercept_m= 1/T * intercept_m Ht[:,:,0] = np.cov(ret[:init,:].T) for j in range(1,T): Ht[:,:,j] = (1-a-b) * intercept_m Ht[:,:,j] = Ht[:,:,j] + a *np.matmul(np.array(ret[j-1,:], ndmin=2).T,np.array(ret[j-1,:], ndmin=2)) Ht[:,:,j] = Ht[:,:,j] + b * Ht[:,:,j-1] return Ht def eq_svech_svd(theta,ret,init,intercept_m, SVD = False): T,N = np.shape(ret) a,b = theta Ht = np.zeros((N, N ,T)) U = np.zeros((N, N ,T)) s = np.zeros((N ,T)) Vh = np.zeros((N, N ,T)) ret_shifted = np.insert(ret, 0, ret[0,:], axis=0) ret_shifted = np.delete(ret_shifted,obj = -1, axis =0) Ht[:,:,0] = np.cov(ret[:init,:].T) if SVD == True: U[:,:,0],s[:,0],Vh[:,:,0] = np.linalg.svd(Ht[:,:,0], full_matrices=True, compute_uv=True, hermitian=True) for n in range(0,N): for m in range(0,N): Ht[n,m,:] = (1-a-b) * intercept_m[n,m]+ a * ret_shifted[:,n]*ret_shifted[:,m] Ht[:,:,0] = np.cov(ret[:init,:].T) for t in range(1,T): Ht[:,:,t] = Ht[:,:,t] + b* Ht[:,:,t-1] if SVD == True: U[:,:,t],s[:,t],Vh[:,:,t] = np.linalg.svd(Ht[:,:,t], full_matrices=True, compute_uv=True, hermitian=True) if SVD==True: return U,s,Vh else: return Ht def loglike_svech(theta,ret,init, intercept_m = None, SVD= False, verbose = False): T,N = np.shape(ret) llf = np.zeros((T,1)) outs =[] if intercept_m is None: intercept_m = np.zeros((N,N)) for t in range(T): intercept_m= intercept_m +np.matmul(np.array(ret[t,:], ndmin=2).T,np.array(ret[t,:], ndmin=2)) intercept_m= 1/T * intercept_m if SVD == True: U,s,Vh = eq_svech_svd(theta,ret,init,intercept_m, SVD) else: Ht = eq_svech_svd(theta,ret,init,intercept_m, SVD) if SVD == True: for i in range(0,T): rV = np.matmul(np.array(ret[i,:], ndmin=2) , Vh[:,:,i].T) Ur = np.matmul(U[:,:,i].T, np.array(ret[i,:], ndmin=2).T) rVS = np.divide(rV, s[:,i]) llf[i] = np.matmul(rVS, Ur) llf = np.sum(llf) + s.sum() else: for i in range(0,T): det = np.log(np.linalg.det(Ht[:,:,i])) mult = np.matmul(np.matmul(np.array(ret[i,:], ndmin=2) , (np.linalg.inv(Ht[:,:,i]))) ,np.array(ret[i,:], ndmin=2).T) llf[i] =(det+ mult) llf = np.sum(llf) if verbose: print(theta,llf) if llf == float("-inf"): #trick to keep alhorithom going even when singular matrix is obtained in some interation llf=0 return llf def composite_pair_loglike_svech(theta,ret,init): T,N = np.shape(ret) Hjt= np.zeros((N, N ,T)) llfj = np.zeros((T,1)) ňintercept_m = np.zeros((N,N)) intercept_m= 1/T*np.matmul(ret.T,ret) Hj = eq_svech_svd(theta,ret,init,intercept_m, False) Hj_inv = np.zeros((N,N,T)) det = np.zeros((T)) det[:] = Hj[0,0,:]* Hj[1,1,:]- Hj[0,1,:]* Hj[1,0,:] Hj_inv[0,0,:] = Hj[1,1,:] /det[:] Hj_inv[1,1,:] = Hj[0,0,:] /det[:] Hj_inv[0,1,:] = -Hj[0,1,:] /det[:] Hj_inv[1,0,:] = -Hj[1,0,:] /det[:] mult = np.einsum('ij, jki, ik ->', ret, Hj_inv, ret) llfj = np.sum(np.log(det))+np.sum(mult) llfj = np.sum(llfj) return llfj def composite_loglike(theta,ret,init, endpair= None, verbose = False): T,N = np.shape(ret) llf = np.zeros((int(N*(N-1)/2),1)) i=0 if endpair is not None: endpair = min(N,endpair) endpair = max(2,endpair) else: endpair = N for m in range(1,endpair): for n in range(0,N): if m+n < N: retj =ret[:,[n,m+n]] llf[i]= composite_pair_loglike_svech(theta,retj,init) i+=1 llf= np.sum(llf) if verbose: print(theta,llf) return llf ----------------------------------------------------------------------------------------------------------- import numpy as np def ore_est_H(a,ret,init): T,N = np.shape(ret) Ht = np.zeros((N, N ,T)) Ht[:,:,0] = np.cov(ret[:init,:].T) for j in range(1,T): Ht[:,:,j] = a * np.exp(-a)*np.matmul(np.array(ret[j-1,:], ndmin=2).T,np.array(ret[j-1,:], ndmin=2)) Ht[:,:,j] = Ht[:,:,j] + np.exp(-a) * Ht[:,:,j-1] return Ht def eq_ore_svd(a,ret,init,SVD= False): T,N = np.shape(ret) Ht = np.zeros((N, N ,T)) U = np.zeros((N, N ,T)) s = np.zeros((N ,T)) Vh = np.zeros((N, N ,T)) ret_shifted = np.insert(ret, 0, ret[0,:], axis=0) ret_shifted = np.delete(ret_shifted,obj = -1, axis =0) Ht[:,:,0] = np.cov(ret[:init,:].T) if SVD== True: U[:,:,0],s[:,0],Vh[:,:,0] = np.linalg.svd(Ht[:,:,0], full_matrices=True, compute_uv=True, hermitian=True) for n in range(0,N): for m in range(0,N): Ht[n,m,:] = a * np.exp(-a)* ret_shifted[:,n]*ret_shifted[:,m] Ht[:,:,0] = np.cov(ret[:init,:].T) for t in range(1,T): Ht[:,:,t] = Ht[:,:,t]+ np.exp(-a) * Ht[:,:,t-1] if SVD == True: U[:,:,t],s[:,t],Vh[:,:,t] = np.linalg.svd(Ht[:,:,t], full_matrices=True, compute_uv=True, hermitian=True) if SVD== True: return U,s,Vh else: return Ht def loglike_ore(a,ret,init, SVD= False, verbose = False): T,N = np.shape(ret) llf = np.zeros((T,1)) outs =[] if SVD== True: U,s,Vh = eq_ore_svd(a,ret,init,SVD) else: Ht = eq_ore_svd(a,ret,init,SVD) if SVD == True: for i in range(0,T): rV = np.matmul(np.array(ret[i,:], ndmin=2) , Vh[:,:,i].T) Ur = np.matmul(U[:,:,i].T, np.array(ret[i,:], ndmin=2).T) rVS = np.divide(rV, s[:,i]) llf[i] = np.matmul(rVS, Ur) #if i % 1000 == 0: #print('llf: {}, a is {}, llf[i] is {}, s is {}'.format(np.sum(llf), a, llf[i], s[:,i].sum())) llf = np.sum(llf) + s.sum() else: for i in range(0,T): det = np.log(np.linalg.det(Ht[:,:,i])) mult = np.matmul(np.matmul(np.array(ret[i,:], ndmin=2) , (np.linalg.inv(Ht[:,:,i]))) ,np.array(ret[i,:], ndmin=2).T) llf[i] =(det+ mult) llf = np.sum(llf) if verbose: print(a,llf) if llf == float("-inf"): llf=0 return llf def composite_pair_loglike_ore(a,ret,init): T,N = np.shape(ret) Hj = np.zeros((N, N ,T)) llfj = np.zeros((T,1)) Hj = eq_ore_svd(a,ret,init, False) Hj_inv = np.zeros((N,N,T)) det = np.zeros((T)) det[:] = Hj[0,0,:]* Hj[1,1,:]- Hj[0,1,:]* Hj[1,0,:] Hj_inv[0,0,:] = Hj[1,1,:] /det[:] Hj_inv[1,1,:] = Hj[0,0,:] /det[:] Hj_inv[0,1,:] = -Hj[0,1,:] /det[:] Hj_inv[1,0,:] = -Hj[1,0,:] /det[:] mult = np.einsum('ij, jki, ik ->', ret, Hj_inv, ret) llfj = np.sum(np.log(det))+np.sum(mult) llfj = np.sum(llfj) return llfj def composite_loglike(a,ret,init, endpair= None,verbose = False): T,N = np.shape(ret) llf = np.zeros((int(N*(N-1)/2),1)) i=0 if endpair is not None: endpair = min(N,endpair) endpair = max(2,endpair) else: endpair = N for m in range(1,endpair): for n in range(0,N): if m+n < N: retj =ret[:,[n,m+n]] llf[i]= composite_pair_loglike_ore(a,retj,init) i+=1 llf= np.sum(llf) if verbose: print(a,llf) return llf ----------------------------------------------------------------------------------------------------------------------------------- { "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import yfinance as yf\n", "import pandas as pd\n", "import os\n", "import pandas_datareader.data as web\n", "import seaborn as sns\n", "import matplotlib.pyplot as plt\n", "from numpy.linalg import cholesky\n", "import random\n", "import math\n", "import numpy as np\n", "import datetime\n", "#%load_ext dotenv\n", "#%dotenv\n", "from IPython.display import display\n", "\n", "import bs4 as bs\n", "import requests\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from arch import arch_model\n", "from scipy.stats import t\n", "from scipy.stats import norm\n", "from scipy.optimize import fmin, minimize,differential_evolution\n", "from scipy import linalg\n", "from scipy.optimize import Bounds\n", "from scipy.ndimage.interpolation import shift\n", "from datetime import datetime\n", "import time" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#all_data_start, all_data_end = '2017-01-01', \"2019-12-31\"\n", "\n", "default_start , default_end = '2010-01-01', \"2021-31-12\"\n", "resp = requests.get('http://en.wikipedia.org/wiki/List_of_S%26P_500_companies')\n", "soup = bs.BeautifulSoup(resp.text, 'lxml')\n", "table = soup.find('table', {'class': 'wikitable sortable'})\n", "tickers = []\n", "date_first_added = []\n", "for row in table.findAll('tr')[1:]:\n", " ticker = row.findAll('td')[0].text\n", " tickers.append(ticker)\n", " date = row.findAll('td')[6].text\n", " date_first_added.append(date)\n", " \n", "\n", "tickers = [s.replace('\\n', '') for s in tickers]\n", "date_first_added = [d.replace('\\n', '') for d in date_first_added]\n", "\n", "date_added = pd.DataFrame(data = (date_first_added), index= tickers)\n", "date_added.columns = [['Date Added']]\n", "date_added = date_added[['Date Added']].apply(lambda x: x[:9], axis = 1)\n", "\n", "#N_ass =len(tickers)\n", "N_ass =500\n", "\n", "start = datetime(2010,1,2)\n", "end = datetime(2022,1,1)\n", "\n", "#start = datetime(2017,1,1)\n", "#end = datetime(2017,1,6)\n", "close_prices_all = yf.download(tickers[:N_ass], start=start, end=end)[['Adj Close','Volume']]\n", "\n", "close_prices = close_prices_all.dropna(axis=1, how='any', thresh=len(close_prices_all)-5, subset=None, inplace=False)\n", "close_prices.shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "N_assets = int(close_prices.shape[1]/2)\n", "mean_volumes = pd.DataFrame(close_prices['Volume'].mean(axis =0))\n", "mean_volumes.columns = ['Mean Volume']\n", "mean_volumes.sort_values(by = 'Mean Volume', ascending = False, inplace =True)\n", "most_traded_indexes = mean_volumes.iloc[:N_assets].index\n", "close_prices_cut = close_prices['Adj Close'][most_traded_indexes]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "#col = close_prices.columns[0]\n", "#a = close_prices[col][close_prices[col].isna()==True].index\n", "close_prices_cut.dropna(axis=0, how='all', thresh=None, subset=None, inplace=True)\n", "if \"APA\" in most_traded_indexes:\n", " close_prices_cut.drop(labels =\"APA\",axis =1, inplace = True)\n", "\n", "if \"GOOGL\" in most_traded_indexes:\n", " close_prices_cut.drop(labels =\"GOOGL\",axis =1, inplace = True) \n", "\n", "Log_ret = pd.DataFrame()\n", "\n", "for col in close_prices_cut.columns:\n", " Log_ret[col] = (np.log(close_prices_cut[col]) - np.log(close_prices_cut[col].shift(1))).dropna(how='all') \n", " \n", "clas_ret = pd.DataFrame()\n", "\n", "for col in close_prices_cut.columns:\n", " clas_ret[col] = ((close_prices_cut[col] / close_prices_cut[col].shift(1))-1).dropna(how='all') " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "close_prices_cut" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Models import and estimation" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#own written modules for candidate models\n", "import DCC_model as DCC\n", "import EWMA_model as EWMA\n", "import ORE_model as ORE\n", "import S_Vech_model as SV # scalar BEKK model\n", "import CCC_model as CCC" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "#Setup parameters\n", "\n", "N_assets = 100\n", "\n", "assets = close_prices_cut.columns[:N_assets]\n", "#df= Log_ret.iloc[:,:N_assets].dropna()\n", "df= clas_ret.iloc[:,:N_assets].dropna()\n", "stocks_names = df.columns\n", "scaling_param = 1\n", "\n", "T,N = np.shape(df)\n", "ret = np.array(df)\n", " \n", "init_step_comp= 5\n", "init_step= 200\n", "\n", "estimate_SVD_versions = True\n", "\n", "Time_cutoff = (500,1000,2000,T) #initial part of the numerical study\n", "\n", "print(T,N,init_step)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#used when the neihbouring composite maximum is deamed the best -\n", "\n", "Time_cutoff =[]\n", "for t in range(2500,T,5):\n", " Time_cutoff.append(t)\n", "if Time_cutoff[-1] != T:\n", " Time_cutoff.append(T)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [], "scrolled": true }, "outputs": [], "source": [ "#Run univariate garch\n", "uni_garch_scale = 100\n", "\n", "model_params, D_diag = DCC.run_univariate_garch(df * uni_garch_scale)\n", "D_diag = D_diag /uni_garch_scale\n", "\n", "D_inv = np.array(pd.DataFrame(D_diag).apply(lambda x: 1/x))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#CCC\n", "R_ccc = CCC.ccceq(ret,D_inv, init_step)\n", "Ht_ccc = np.zeros((N,N,T))\n", "for t in range(T):\n", " Ht_ccc[:,:,t] = np.matmul(np.matmul(np.diagflat(D_diag[t,:]),R_ccc[:,:,t]),np.diagflat(D_diag[t,:]))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [], "scrolled": true }, "outputs": [], "source": [ "DCC_svd_est=[]\n", "for cut in Time_cutoff:#DCC SVD\n", " if estimate_SVD_versions == True:\n", " # lambda functions definds that the following ineq must hold \n", " cons = ({'type': 'ineq', 'fun': lambda x: -x[0] -x[1] +1})\n", " bnds = ((0, 0.04), (0.95, 0.9997))\n", "\n", " %time opt_out_DCC_SVD = minimize(DCC.loglike_dcc_svd, np.array([0.001, 0.98]), args = (ret[:cut],D_inv[:cut],init_step,True,True),\\\n", " method='Nelder-Mead', bounds=bnds,tol = (0.1))\n", " print(\"DCC SVD estimates alpha: {}, beta: {}\".format(opt_out_DCC_SVD.x[0],opt_out_DCC_SVD.x[1]))\n", " DCC_svd_est.append(opt_out_DCC_SVD.x)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [], "scrolled": true }, "outputs": [], "source": [ "DCC_est=[]\n", "for cut in Time_cutoff:#DCC\n", " # lambda functions definds that the following ineq must hold \n", " cons = ({'type': 'ineq', 'fun': lambda x: -x[0] -x[1] +1})\n", " bnds = ((0, 0.1), (0.6, 0.9997))\n", "\n", " %time opt_out_DCC = minimize(DCC.loglike_dcc_svd, np.array([0.001, 0.99]), args = (ret[:cut],D_inv[:cut],init_step,False,True),\\\n", " method='nelder-mead', bounds=bnds, constraints=cons)\n", " print(\"DCC SVD estimates alpha: {}, beta: {}\".format(opt_out_DCC_SVD.x[0],opt_out_DCC_SVD.x[1]))\n", " DCC_est.append(opt_out_DCC.x)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "#DCC composite ML\n", "# lambda functions definds that the following ineq must hold \n", "DCC_comp_est=[np.array([0.04, 0.88])]\n", "for cut in Time_cutoff:\n", " cons = ({'type': 'ineq', 'fun': lambda x: -x[0] -x[1] +1})\n", " bnds = ((0, 0.1), (0.6, 0.9997))\n", "\n", " %time opt_out_composite_DCC = minimize(DCC.composite_loglike,DCC_comp_est[-1], method='nelder-mead', \\\n", " args = (ret[:cut],D_inv[:cut],init_step_comp,2), bounds=bnds, constraints=cons)\n", " print(\"DCC composite ML estimates alpha: {}, beta: {}\".format(opt_out_composite_DCC.x[0],opt_out_composite_DCC.x[1]))\n", " DCC_comp_est.append(opt_out_composite_DCC.x)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "#DCC composite ML full\n", "# lambda functions definds that the following ineq must hold \n", "DCC_comp_est_full=[] \n", "for cut in Time_cutoff:\n", " cons = ({'type': 'ineq', 'fun': lambda x: -x[0] -x[1] +1})\n", " bnds = ((0, 0.1), (0.6, 0.9997))\n", "\n", " %time opt_out_composite_DCC = minimize(DCC.composite_loglike, np.array([0.04, 0.88]), method='nelder-mead', \\\n", " args = (ret[:cut],D_inv[:cut],init_step_comp,N), bounds=bnds, constraints=cons)\n", " print(\"DCC composite ML estimates alpha: {}, beta: {}\".format(opt_out_composite_DCC.x[0],opt_out_composite_DCC.x[1]))\n", " DCC_comp_est_full.append(opt_out_composite_DCC.x)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [], "scrolled": true }, "outputs": [], "source": [ "#EWMA SVD\n", "EWMA_svd_est=[]\n", "for cut in Time_cutoff:\n", " if estimate_SVD_versions == True:\n", " # lambda functions definds that the following ineq must hold \n", " bnds = ((0,0.1),)\n", "\n", " %time opt_out_ewma_svd = minimize(EWMA.loglike_ewma_svd,x0=0.0015, bounds=bnds, args = (ret[:cut],init_step,True))\n", " print(\"EWMA SVD estimates alpha: {}\".format(opt_out_ewma_svd.x))\n", " EWMA_svd_est.append(opt_out_ewma_svd.x)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [], "scrolled": true }, "outputs": [], "source": [ "#EWMA \n", "bnds = ((0,0.1),)\n", "EWMA_est=[]\n", "for cut in Time_cutoff:\n", " %time opt_out_ewma = minimize(EWMA.loglike_ewma_svd,x0=0.0015, bounds=bnds, method='nelder-mead', args = (ret[:cut],init_step,False,True))\n", " print(\"EWMA estimates alpha: {}\".format(opt_out_ewma.x))\n", " EWMA_est.append(opt_out_ewma.x)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "#EWMA composite ML\n", "bnds = ((0,0.1),)\n", "EWMA_comp_est=[np.array([0.025])]\n", "for cut in Time_cutoff:\n", " %time opt_out_composite_ewma = minimize(EWMA.composite_loglike,x0=EWMA_comp_est[-1], method='nelder-mead',bounds=bnds, args = (ret[:cut],init_step_comp,2))\n", " print(\"EWMA composite ML estimates alpha: {}\".format(opt_out_composite_ewma.x))\n", " EWMA_comp_est.append(opt_out_composite_ewma.x)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#EWMA composite ML full\n", "bnds = ((0,0.1),)\n", "EWMA_comp_est_full=[]\n", "for cut in Time_cutoff:\n", " %time opt_out_composite_ewma = minimize(EWMA.composite_loglike,x0=0.025, method='nelder-mead',bounds=bnds, args = (ret[:cut],init_step_comp,N))\n", " print(\"EWMA composite ML estimates alpha: {}\".format(opt_out_composite_ewma.x))\n", " EWMA_comp_est_full.append(opt_out_composite_ewma.x)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [], "scrolled": true }, "outputs": [], "source": [ "#ORE SVD\n", "ORE_svd_est=[]\n", "for cut in Time_cutoff:\n", " if estimate_SVD_versions == True:\n", " # lambda functions definds that the following ineq must hold \n", " bnds = ((0,0.1),)\n", "\n", " %time opt_out_ore_svd = minimize(ORE.loglike_ore,x0=0.04,bounds=bnds, args = (ret[:cut],init_step,True))\n", "\n", "\n", " print(\"ORE SVD estimates alpha: {}\".format(opt_out_ore_svd.x))\n", " ORE_svd_est.append(opt_out_ore_svd.x)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "#ORE\n", "bnds = ((0,0.1),)\n", "ORE_est=[]\n", "for cut in Time_cutoff:\n", "\n", " %time opt_out_ore = minimize(ORE.loglike_ore,x0=0.04,method='nelder-mead', bounds=bnds, args = (ret[:cut],init_step,False))\n", " print(\"ORE estimates alpha: {}\".format(opt_out_ore.x))\n", " ORE_est.append(opt_out_ore.x)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "#ORE composite ML\n", "bnds = ((0,0.1),)\n", "ORE_comp_est=[np.array(0.025)]\n", "for cut in Time_cutoff:\n", " %time opt_out_ore_composite = minimize(ORE.composite_loglike,x0=0.025, bounds=bnds,method='nelder-mead', args = (ret[:cut],init_step_comp,2))\n", " print(\"ORE composite ML estimates alpha: {}\".format(opt_out_ore_composite.x))\n", " ORE_comp_est.append(opt_out_ore_composite.x)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#ORE composite ML full\n", "bnds = ((0,0.1),)\n", "ORE_comp_est_full=[]\n", "for cut in Time_cutoff:\n", " %time opt_out_ore_composite = minimize(ORE.composite_loglike,x0=0.025, bounds=bnds,method='nelder-mead', args = (ret[:cut],init_step_comp,N))\n", " print(\"ORE composite ML estimates alpha: {}\".format(opt_out_ore_composite.x))\n", " ORE_comp_est_full.append(opt_out_ore_composite.x)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#SVECH SVD\n", "SVECH_svd_est=[]\n", "for cut in Time_cutoff[2:]:\n", " if estimate_SVD_versions == True:\n", " cons = ({'type': 'ineq', 'fun': lambda x: -x[0] -x[1] +1})\n", " bnds = ((0, 0.1), (0.6, 0.99))\n", "\n", " %time opt_out_svech_svd = minimize(SV.loglike_svech,x0=np.array([0.0009, 0.99]), bounds = bnds, constraints =cons, tol=0.0001, args = (ret[:cut],init_step, None, True,True))\n", " print(\"SVECH SVD estimates alpha: {}, beta: {}\".format(opt_out_svech_svd.x[0],opt_out_svech_svd.x[1]))\n", " SVECH_svd_est.append(opt_out_svech_svd.x)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "#SVECH SVD\n", "SVECH_svd_est=[]\n", "for cut in Time_cutoff:\n", " if estimate_SVD_versions == True:\n", " cons = ({'type': 'ineq', 'fun': lambda x: -x[0] -x[1] +1})\n", " bnds = ((0, 0.1), (0.6, 0.99))\n", "\n", " %time opt_out_svech_svd = minimize(SV.loglike_svech,x0=np.array([0.004, 0.9699]), bounds = bnds, constraints =cons, tol=0.0001, args = (ret[:cut],init_step, None, True,True))\n", " print(\"SVECH SVD estimates alpha: {}, beta: {}\".format(opt_out_svech_svd.x[0],opt_out_svech_svd.x[1]))\n", " SVECH_svd_est.append(opt_out_svech_svd.x)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "#SVECH\n", "SVECH_est=[]\n", "for cut in Time_cutoff:\n", " cons = ({'type': 'ineq', 'fun': lambda x: -x[0] -x[1] +1})\n", " bnds = ((0, 0.1), (0, 0.99))\n", "\n", " %time opt_out_svech = minimize(SV.loglike_svech,x0=np.array([0.002, 0.96]), bounds=bnds,method='nelder-mead', constraints=cons, args = (ret[:cut],init_step, None, False ))\n", " print(\"SVECH estimates alpha: {}, beta: {}\".format(opt_out_svech.x[0],opt_out_svech.x[1]))\n", " SVECH_est.append(opt_out_svech.x)\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "#SVECH composite ML\n", "\n", "SVECH_comp_est=[np.array([0.06, 0.90])]\n", "for cut in Time_cutoff:\n", " %time opt_out_composite_svech = minimize(SV.composite_loglike,x0=SVECH_comp_est[-1],method='nelder-mead', args = (ret[:cut],init_step_comp,2))\n", " print(\"SVECH compososite ML estimates alpha: {}, beta: {}\".format(opt_out_composite_svech.x[0],opt_out_composite_svech.x[1]))\n", " SVECH_comp_est.append(opt_out_composite_svech.x)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#SVECH composite ML full\n", "\n", "SVECH_comp_est_full=[]\n", "for cut in Time_cutoff:\n", " %time opt_out_composite_svech = minimize(SV.composite_loglike,x0=np.array([0.06, 0.90]),method='nelder-mead',bounds=bnds ,constraints=cons, args = (ret[:cut],init_step_comp,N))\n", " print(\"SVECH compososite ML estimates alpha: {}, beta: {}\".format(opt_out_composite_svech.x[0],opt_out_composite_svech.x[1]))\n", " SVECH_comp_est_full.append(opt_out_composite_svech.x)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#comment out the models which parameters were not estimated\n", "\n", "model_params =[\n", "#EWMA_est,\n", "#EWMA_svd_est,\n", "EWMA_comp_est,\n", "#EWMA_comp_est_full,\n", "#DCC_est,\n", "#DCC_svd_est,\n", "DCC_comp_est,\n", "#DCC_comp_est_full,\n", "#SVECH_est,\n", "#SVECH_svd_est,\n", "SVECH_comp_est,\n", "#SVECH_comp_est_full,\n", "#ORE_est,\n", "#ORE_svd_est,\n", "ORE_comp_est,\n", "#ORE_comp_est_full\n", "]\n", "estimates_names =[\n", "#'EWMA_est',\n", "#'EWMA_svd_est',\n", "'EWMA_comp_est',\n", "#'EWMA_comp_est_full',\n", "#'DCC_est',\n", "#'DCC_svd_est',\n", "'DCC_comp_est',\n", "#'DCC_comp_est_full',\n", "#'SVECH_est',\n", "#'SVECH_svd_est',\n", "'SVECH_comp_est',\n", "#'SVECH_comp_est_full',\n", "#'ORE_est',\n", "#'ORE_svd_est',\n", "'ORE_comp_est',\n", "#'ORE_comp_est_full'\n", "]\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "estimates = pd.DataFrame(columns=Time_cutoff)\n", "\n", "\n", "for model in model_params:\n", " temp = pd.DataFrame(columns=Time_cutoff)\n", " for i,col in enumerate(estimates.columns):\n", " temp[col] =[model[i]]\n", " \n", " estimates = pd.concat([estimates,temp])\n", "\n", "\n", "estimates['ind'] = estimates_names\n", "estimates.set_index('ind')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "estimates['ind'] = estimates_names\n", "estimates.set_index('ind')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "estimates.to_excel('estimates.xlsx')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "estimates.to_csv(\"csvestimates.csv\", sep='&')" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": { "height": "calc(100% - 180px)", "left": "10px", "top": "150px", "width": "299.2px" }, "toc_section_display": true, "toc_window_display": true }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 4 } -------------------------------------------------------------------------------------------------------------- { "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import yfinance as yf\n", "import pandas as pd\n", "import os\n", "import pandas_datareader.data as web\n", "import seaborn as sns\n", "import matplotlib.pyplot as plt\n", "from numpy.linalg import cholesky\n", "import random\n", "import math\n", "import numpy as np\n", "import datetime\n", "#%load_ext dotenv\n", "#%dotenv\n", "from IPython.display import display\n", "\n", "import bs4 as bs\n", "import requests\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from arch import arch_model\n", "from scipy.stats import t\n", "from scipy.stats import norm\n", "from scipy.optimize import fmin, minimize,differential_evolution\n", "from scipy import linalg\n", "from scipy.optimize import Bounds\n", "from scipy.ndimage.interpolation import shift\n", "from datetime import datetime\n", "import time\n", "import tqdm\n", "import itertools" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#all_data_start, all_data_end = '2017-01-01', \"2019-12-31\"\n", "\n", "default_start , default_end = '2010-01-01', \"2021-31-12\"\n", "resp = requests.get('http://en.wikipedia.org/wiki/List_of_S%26P_500_companies')\n", "soup = bs.BeautifulSoup(resp.text, 'lxml')\n", "table = soup.find('table', {'class': 'wikitable sortable'})\n", "tickers = []\n", "date_first_added = []\n", "for row in table.findAll('tr')[1:]:\n", " ticker = row.findAll('td')[0].text\n", " tickers.append(ticker)\n", " date = row.findAll('td')[6].text\n", " date_first_added.append(date)\n", " \n", "\n", "tickers = [s.replace('\\n', '') for s in tickers]\n", "date_first_added = [d.replace('\\n', '') for d in date_first_added]\n", "\n", "date_added = pd.DataFrame(data = (date_first_added), index= tickers)\n", "date_added.columns = [['Date Added']]\n", "date_added = date_added[['Date Added']].apply(lambda x: x[:9], axis = 1)\n", "\n", "N_ass =len(tickers)\n", "#N_ass =500\n", "\n", "start = datetime(2010,1,2)\n", "end = datetime(2022,1,1)\n", "\n", "#start = datetime(2017,1,1)\n", "#end = datetime(2017,1,6)\n", "close_prices_all = yf.download(tickers[:N_ass], start=start, end=end)[['Adj Close','Volume']]\n", "SP500 = yf.download('^GSPC', start=start, end=end)[['Adj Close']]\n", "\n", "close_prices = close_prices_all.dropna(axis=1, how='any', thresh=len(close_prices_all)-5, subset=None, inplace=False)\n", "close_prices.shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "N_assets = int(close_prices.shape[1]/2)\n", "mean_volumes = pd.DataFrame(close_prices['Volume'].mean(axis =0))\n", "mean_volumes.columns = ['Mean Volume']\n", "mean_volumes.sort_values(by = 'Mean Volume', ascending = False, inplace =True)\n", "most_traded_indexes = mean_volumes.iloc[:N_assets].index\n", "close_prices_cut = close_prices['Adj Close'][most_traded_indexes]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "#col = close_prices.columns[0]\n", "#a = close_prices[col][close_prices[col].isna()==True].index\n", "close_prices_cut.dropna(axis=0, how='all', thresh=None, subset=None, inplace=True)\n", "if \"APA\" in most_traded_indexes:\n", " close_prices_cut.drop(labels =\"APA\",axis =1, inplace = True)\n", "\n", "if \"GOOGL\" in most_traded_indexes:\n", " close_prices_cut.drop(labels =\"GOOGL\",axis =1, inplace = True) \n", " \n", "Log_ret = pd.DataFrame()\n", "\n", "for col in close_prices_cut.columns:\n", " Log_ret[col] = (np.log(close_prices_cut[col]) - np.log(close_prices_cut[col].shift(1))).dropna(how='all') \n", " \n", "clas_ret = pd.DataFrame()\n", "\n", "for col in close_prices_cut.columns:\n", " clas_ret[col] = ((close_prices_cut[col] / close_prices_cut[col].shift(1))-1).dropna(how='all') \n", " \n", "SP_ret = np.array(((SP500 / SP500.shift(1))-1).dropna(how='all'))\n", "dates= clas_ret.index" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "\n", "clas_ret" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "corr = close_prices_cut.corr()\n", "(np.isnan(corr[(abs(corr)>0.99) &(abs(corr)<1)]).sum()-close_prices_cut.shape[1]).sum()\n", "#f = plt.figure(figsize=(19, 15))\n", "#plt.matshow(df.corr(), fignum=f.number)\n", "#plt.xticks(range(df.select_dtypes(['number']).shape[1]), df.select_dtypes(['number']).columns, fontsize=14, rotation=45)\n", "#plt.yticks(range(df.select_dtypes(['number']).shape[1]), df.select_dtypes(['number']).columns, fontsize=14)\n", "#cb = plt.colorbar()\n", "#cb.ax.tick_params(labelsize=14)\n", "#plt.title('Correlation Matrix', fontsize=16);" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "corr[(abs(corr)>0.995) &(abs(corr)<1)].dropna(axis =1, how ='all').dropna(axis =0, how ='all')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Models import and estimation\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#own written modules for candidate models\n", "import DCC_model as DCC\n", "import EWMA_model as EWMA\n", "import ORE_model as ORE\n", "import S_Vech_model as SV # scalar BEKK model\n", "import CCC_model as CCC" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "#Setup parameters\n", "\n", "N_assets = 100\n", "\n", "assets = close_prices_cut.columns[:N_assets]\n", "#df= Log_ret.iloc[:,:N_assets].dropna()\n", "df= clas_ret.iloc[:,:N_assets].dropna()\n", "stocks_names = df.columns\n", "scaling_param = 1\n", "\n", "T,N = np.shape(df)\n", "ret = np.array(df)\n", " \n", "valid_wind = 300\n", "init_step=200\n", "init_step_comp=5\n", "\n", "estimate_SVD_versions = True\n", "\n", "\n", "\n", "print(T,N,init_step)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [], "scrolled": true }, "outputs": [], "source": [ "#Run univariate garch\n", "uni_garch_scale = 100\n", "\n", "model_params, D_diag = DCC.run_univariate_garch(df * uni_garch_scale)\n", "D_diag = D_diag /uni_garch_scale\n", "\n", "D_inv = np.array(pd.DataFrame(D_diag).apply(lambda x: 1/x))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Estimating the parameters of the modes - Not needed when a file with parameter estimates is available\n", "Used before estimating the models for different time period when the thesis was just being implemented. If the parameters are estimated usign Thesis_code-set_up this part of the code can be skipped to $ \\textbf{2. Construct the matrices from excel files}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#CCC\n", "R_ccc = CCC.ccceq(ret,D_inv,200)\n", "Ht_ccc = np.zeros((N,N,T))\n", "for t in range(T):\n", " Ht_ccc[:,:,t] = np.matmul(np.matmul(np.diagflat(D_diag[t,:]),R_ccc[:,:,t]),np.diagflat(D_diag[t,:]))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [], "scrolled": true }, "outputs": [], "source": [ "#DCC SVD\n", "if estimate_SVD_versions == True:\n", "# lambda functions definds that the following ineq must hold \n", " cons = ({'type': 'ineq', 'fun': lambda x: -x[0] -x[1] +1})\n", " bnds = ((0, 0.1), (0.6, 0.9997))\n", "\n", " %time opt_out_DCC_SVD = minimize(DCC.loglike_dcc_svd, np.array([0.00028, 0.999]), args = (ret,D_inv,init_step,True,True),\\\n", " method='nelder-mead', bounds=bnds, constraints=cons)\n", " print(\"DCC SVD estimates alpha: {}, beta: {}\".format(opt_out_DCC_SVD.x[0],opt_out_DCC_SVD.x[1]))\n", "\n", " R= DCC.dcceq(opt_out_DCC_SVD.x, ret, D_inv,init_step)\n", " Ht_dcc_svd = np.zeros((N,N,T))\n", "\n", " for t in range(T):\n", " Ht_dcc_svd[:,:,t] = np.matmul(np.matmul(np.diagflat(D_diag[t,:]),R[:,:,t]),np.diagflat(D_diag[t,:]))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [], "scrolled": true }, "outputs": [], "source": [ "#DCC\n", "# lambda functions definds that the following ineq must hold \n", "cons = ({'type': 'ineq', 'fun': lambda x: -x[0] -x[1] +1})\n", "bnds = ((0, 0.1), (0.6, 0.9997))\n", "\n", "%time opt_out_DCC = minimize(DCC.loglike_dcc_svd, np.array([0.01, 0.95]), args = (ret,D_inv,init_step,False),\\\n", " method='nelder-mead', bounds=bnds, constraints=cons)\n", "print(\"DCC estimates alpha: {}, beta: {}\".format(opt_out_DCC.x[0],opt_out_DCC.x[1]))\n", "\n", "R= DCC.dcceq(opt_out_DCC.x, ret,D_inv,init_step)\n", "Ht_dcc = np.zeros((N,N,T))\n", "for t in range(T):\n", " Ht_dcc[:,:,t] = np.matmul(np.matmul(np.diagflat(D_diag[t,:]),R[:,:,t]),np.diagflat(D_diag[t,:]))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "#DCC composite ML\n", "# lambda functions definds that the following ineq must hold \n", "cons = ({'type': 'ineq', 'fun': lambda x: -x[0] -x[1] +1})\n", "bnds = ((0, 0.1), (0.6, 0.9997))\n", "\n", "%time opt_out_composite_DCC = minimize(DCC.composite_loglike, np.array([0.04, 0.88]), method='nelder-mead', \\\n", " args = (ret,D_inv,init_step_comp,2), bounds=bnds, constraints=cons)\n", "print(\"DCC composite ML estimates alpha: {}, beta: {}\".format(opt_out_composite_DCC.x[0],opt_out_composite_DCC.x[1]))\n", "\n", "R_comp= DCC.dcceq(opt_out_composite_DCC.x, ret, D_inv,init_step)\n", "Ht_comp_dcc = np.zeros((N,N,T))\n", "for t in range(T):\n", " Ht_comp_dcc[:,:,t] = np.matmul(np.matmul(np.diagflat(D_diag[t,:]),R_comp[:,:,t]),np.diagflat(D_diag[t,:]))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [], "scrolled": true }, "outputs": [], "source": [ "#EWMA SVD\n", "if estimate_SVD_versions == True:\n", " # lambda functions definds that the following ineq must hold \n", " bnds = ((0,0.1),)\n", "\n", " %time opt_out_ewma_svd = minimize(EWMA.loglike_ewma_svd,x0=0.05, bounds=bnds, args = (ret,init_step,True))\n", " print(\"EWMA SVD estimates alpha: {}\".format(opt_out_ewma_svd.x))\n", "\n", " Ht_ewma_svd = np.zeros((N,N,T))\n", " Ht_ewma_svd = EWMA.ewma_est_H(opt_out_ewma_svd.x,ret,init_step)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [], "scrolled": true }, "outputs": [], "source": [ "#EWMA \n", "bnds = ((0,0.1),)\n", "\n", "%time opt_out_ewma = minimize(EWMA.loglike_ewma_svd,x0=0.001, bounds=bnds, args = (ret,init_step,False,True))\n", "print(\"EWMA estimates alpha: {}\".format(opt_out_ewma.x))\n", "\n", "Ht_ewma= np.zeros((N,N,T))\n", "Ht_ewma = EWMA.ewma_est_H(opt_out_ewma.x,ret,init_step)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#EWMA composite ML\n", "bnds = ((0,0.1),)\n", "\n", "%time opt_out_composite_ewma = minimize(EWMA.composite_loglike,x0=0.025, method='nelder-mead',bounds=bnds, args = (ret,init_step_comp,2))\n", "print(\"EWMA composite ML estimates alpha: {}\".format(opt_out_composite_ewma.x))\n", "\n", "Ht_comp_ewma= np.zeros((N,N,T))\n", "Ht_comp_ewma = EWMA.ewma_est_H(opt_out_composite_ewma.x,ret,init_step)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [], "scrolled": true }, "outputs": [], "source": [ "#ORE SVD\n", "if estimate_SVD_versions == True:\n", " # lambda functions definds that the following ineq must hold \n", " bnds = ((0,0.1),)\n", "\n", " %time opt_out_ore_svd = minimize(ORE.loglike_ore,x0=0.04, bounds=bnds, args = (ret,init_step,True))\n", "\n", "\n", " print(\"ORE SVD estimates alpha: {}\".format(opt_out_ore_svd.x))\n", "\n", " Ht_ore_svd = np.zeros((N,N,T))\n", " Ht_ore_svd = ORE.ore_est_H(opt_out_ore_svd.x,ret,init_step)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [], "scrolled": true }, "outputs": [], "source": [ "#ORE\n", "bnds = ((0,0.02),)\n", "\n", "%time opt_out_ore = minimize(ORE.loglike_ore,x0=0.02,method ='nelder-mead', bounds=bnds, args = (ret,init_step,False,True))\n", "\n", "\n", "print(\"ORE estimates alpha: {}\".format(opt_out_ore.x))\n", "\n", "Ht_ore = np.zeros((N,N,T))\n", "Ht_ore = ORE.ore_est_H(opt_out_ore.x,ret,init_step)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "#ORE composite ML\n", "bnds = ((0,0.1),)\n", "\n", "%time opt_out_ore_composite = minimize(ORE.composite_loglike,x0=0.025, bounds=bnds,method='nelder-mead', args = (ret,init_step_comp,2, True))\n", "\n", "\n", "print(\"ORE composite ML estimates alpha: {}\".format(opt_out_ore_composite.x))\n", "\n", "Ht_comp_ore = np.zeros((N,N,T))\n", "Ht_comp_ore = ORE.ore_est_H(opt_out_ore_composite.x,ret,init_step)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#ORE composite ML all pairs\n", "bnds = ((0,0.1),)\n", "\n", "%time opt_out_ore_composite = minimize(ORE.composite_loglike,x0=0.025, bounds=bnds,method='nelder-mead', args = (ret,init_step_comp,N))\n", "\n", "\n", "print(\"ORE composite ML estimates alpha: {}\".format(opt_out_ore_composite.x))\n", "\n", "Ht_comp_ore = np.zeros((N,N,T))\n", "Ht_comp_ore = ORE.ore_est_H(opt_out_ore_composite.x,ret,init_step)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [], "scrolled": true }, "outputs": [], "source": [ "#SVECH SVD\n", "if estimate_SVD_versions == True:\n", " intercept_m= np.zeros((N,N))\n", " for t in range(T):\n", " intercept_m= intercept_m +np.matmul(np.array(ret[t,:], ndmin=2).T,np.array(ret[t,:], ndmin=2))\n", " intercept_m= 1/T * intercept_m \n", "\n", " cons = ({'type': 'ineq', 'fun': lambda x: -x[0] -x[1] +1})\n", " bnds = ((0, 0.1), (0.6, 0.99))\n", "\n", " %time opt_out_svech_svd = minimize(SV.loglike_svech,x0=np.array([0.002, 0.97]),method='nelder-mead', bounds=bnds, constraints=cons, args = (ret,init_step,intercept_m, True,True))\n", " print(\"SVECH SVD estimates alpha: {}, beta: {}\".format(opt_out_svech_svd.x[0],opt_out_svech_svd.x[1]))\n", "\n", " Ht_svech_svd = np.zeros((N,N,T))\n", " Ht_svech_svd = svech_est_H(SV.opt_out_svech_svd.x,ret,init_step,intercept_m)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "#SVECH\n", "intercept_m= np.zeros((N,N))\n", "\n", "for t in range(T):\n", " intercept_m= intercept_m +np.matmul(np.array(ret[t,:], ndmin=2).T,np.array(ret[t,:], ndmin=2))\n", "intercept_m= 1/T * intercept_m \n", "\n", "cons = ({'type': 'ineq', 'fun': lambda x: -x[0] -x[1] +1})\n", "bnds = ((0, 0.1), (0, 0.99))\n", "\n", "%time opt_out_svech = minimize(SV.loglike_svech,x0=np.array([0.002333, 0.96825828]),method='nelder-mead', bounds=bnds, constraints=cons, args = (ret[:,:],300,None, True,True))\n", "print(\"SVECH estimates alpha: {}, beta: {}\".format(opt_out_svech.x[0],opt_out_svech.x[1]))\n", "\n", "Ht_svech= np.zeros((N,N,T))\n", "Ht_svech = SV.svech_est_H(opt_out_svech.x,ret,init_step,intercept_m)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "#SVECH composite ML\n", "cons = ({'type': 'ineq', 'fun': lambda x: -x[0] -x[1] +1})\n", "bnds = ((0, 0.1), (0.5, 0.99))\n", "\n", "%time opt_out_composite_svech = minimize(SV.composite_loglike,x0=np.array([0.06, 0.90]),method='nelder-mead',bounds=bnds ,constraints=cons, args = (ret,init_step_comp,2,True))\n", "print(\"SVECH compososite ML estimates alpha: {}, beta: {}\".format(opt_out_composite_svech.x[0],opt_out_composite_svech.x[1]))\n", "\n", "intercept_m = np.zeros((N,N))\n", "for t in range(T):\n", " intercept_m= intercept_m +np.matmul(np.array(ret[t,:], ndmin=2).T,np.array(ret[t,:], ndmin=2))\n", "intercept_m= 1/T * intercept_m \n", "\n", "Ht_comp_svech= np.zeros((N,N,T))\n", "Ht_comp_svech = SV.svech_est_H(opt_out_composite_svech.x,ret,init_step,intercept_m)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Construct the matrices from excel files" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "params = pd.read_excel('N100v2.xlsx',index_col='ind')\n", "if 'Unnamed: 0' in params.columns:\n", " params.drop(labels = 'Unnamed: 0', inplace = True, axis = 1)\n", "params" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "params = pd.read_excel('text_outputs.xlsx',index_col='ind',sheet_name='output')\n", "if 'Unnamed: 0' in params.columns:\n", " params.drop(labels = 'Unnamed: 0', inplace = True, axis = 1)\n", "params" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ests = []\n", "for row in params.iterrows():\n", " est =row[1]\n", " est = [e.replace(\"[\",\"\") for e in est]\n", " est = [e.replace(\"]\",\"\") for e in est]\n", " \n", " ests.append(est)\n", "b= pd.DataFrame(ests)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for m in range(params.shape[0]):\n", " for n in range(params.shape[1]):\n", " c =list(filter(len, b.iloc[m,n].split(\" \")))\n", " c = list(map(float, c)) \n", " c = np.array(list(map(np.array, c)))\n", " params.iloc[m,n]= c\n", "params" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Time_cutoff =[]\n", "if params.shape[1]==4:\n", " Time_cutoff = (500,1000,2000,T)\n", "else: \n", " for t in range(300,T,5):\n", " Time_cutoff.append(t)\n", " if Time_cutoff[-1] != T:\n", " Time_cutoff.append(T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## chart sizesettings " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ " SMALL_SIZE = 8\n", " MEDIUM_SIZE = 12\n", " BIGGER_SIZE = 25 # BIGGER_SIZE = 25\n", " \n", " \n", " plt.rc('font', size=BIGGER_SIZE) # controls default text sizes\n", " plt.rc('axes', titlesize=BIGGER_SIZE) # fontsize of the axes title\n", " plt.rc('axes', labelsize=BIGGER_SIZE) # fontsize of the x and y labels\n", " plt.rc('xtick', labelsize=BIGGER_SIZE) # fontsize of the tick labels\n", " plt.rc('ytick', labelsize=BIGGER_SIZE) # fontsize of the tick labels\n", " plt.rc('legend', fontsize=MEDIUM_SIZE) # legend fontsize\n", " plt.rc('figure', titlesize=40) # fontsize of the figure title" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "#Create figure used in thesis, might not work\n", " Models = ['EWMA',\"ORE\", 'BEKK', 'DCC']\n", " #method =['ML','SVD', 'NCLM','FCLM']\n", " method =['NCLM']\n", " line_s=['solid','dotted','dashed','dashdot']\n", " line_c=['b','m','g','k']#['purple','blue','skyblue','deepskyblue']\n", "\n", " mx= 0.07\n", "\n", " fig, axes = plt.subplots(nrows=2, ncols=2,figsize=(25,20))\n", " for m in range(2):\n", " for n in range (2):\n", " axes[m,n].set_title(Models[n+2*m])\n", " if m*n!=1:\n", " axes[m,n].set_ylabel(r\"$\\alpha$\" , color=\"black\")\n", " for i in range(1):\n", " x=dates[range(295,T,5)]\n", " #y = params.iloc[(m+2*n)*4+i]\n", " y= params.iloc[(m*2+n)]\n", " # for p in range(len(y[0])):\n", " # for p in range(len(y[0,0])):\n", " for p in range(len(y.iloc[0])):\n", " if p==0: \n", " par = np.array([item[p] for item in y]) \n", " a = axes[m,n].plot(x,par,color=line_c[i], lw=5, ls='-', label=(r\"$\\alpha$ \"))\n", " axes[m,n].set_ylim(0,mx*1.5) \n", " if p ==1:\n", " par = np.array([item[p] for item in y]) \n", " if i ==0:\n", " ax2= axes[m,n].twinx()\n", " if n==1:\n", " ax2.set_ylabel(r\"$\\beta$\" , color=\"black\")\n", " ax2.set_ylim(par.min()*0.8,par.max()*1.2) \n", " b= ax2.plot(x,par,color=line_c[i+1], lw=5, ls='--', label=(r\"$\\beta$ \"))\n", " lns=a+b \n", " labs = [l.get_label() for l in lns]\n", "\n", " if p==0:\n", " axes[m,n].legend(loc = 'upper right')\n", " else:\n", " axes[m,n].legend(lns, labs, loc = 'upper right') \n", " #ax2.legend(loc = 'lower right') \n", " axes[m,n].set_xlabel(\"year\" , color=\"black\")\n", "\n", " #fig.savefig('paramest.pdf')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "#setup matrices, comment out unwanted matrices\n", "Ht_ewma = np.zeros((N,N,T)) \n", "Ht_ewma_svd = np.zeros((N,N,T))\n", "Ht_comp_ewma = np.zeros((N,N,T)) \n", "Ht_comp_ewma_full = np.zeros((N,N,T)) \n", "\n", "Rt_dcc= np.zeros((N,N,T)) \n", "Rt_dcc_svd = np.zeros((N,N,T))\n", "Rt_comp_dcc= np.zeros((N,N,T)) \n", "Rt_comp_dcc_full= np.zeros((N,N,T)) \n", "Ht_dcc= np.zeros((N,N,T)) \n", "Ht_dcc_svd = np.zeros((N,N,T))\n", "Ht_comp_dcc= np.zeros((N,N,T)) \n", "Ht_comp_dcc_full= np.zeros((N,N,T)) \n", "\n", "Ht_svech= np.zeros((N,N,T)) \n", "Ht_svech_svd = np.zeros((N,N,T))\n", "Ht_comp_svech= np.zeros((N,N,T)) \n", "Ht_comp_svech_full= np.zeros((N,N,T)) \n", "\n", "Ht_ore= np.zeros((N,N,T)) \n", "Ht_ore_svd = np.zeros((N,N,T))\n", "Ht_comp_ore= np.zeros((N,N,T)) \n", "Ht_comp_ore_full= np.zeros((N,N,T)) \n", "\n", "estimates_names =[\n", "'EWMA_est',\n", "'EWMA_svd_est',\n", "'EWMA_comp_est',\n", "'DCC_est',\n", "'DCC_svd_est',\n", "'DCC_comp_est',\n", "'SVECH_est',\n", "'SVECH_svd_est',\n", "'SVECH_comp_est',\n", "'ORE_est',\n", "'ORE_svd_est',\n", "'ORE_comp_est'\n", "]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for a,t in enumerate(Time_cutoff):\n", " if t == Time_cutoff[0]:\n", " s=0\n", " else:\n", " s= Time_cutoff[0]\n", " \n", " print(i,s,t) \n", " Ht_ewma[:,:,s:t] = EWMA.ewma_est_H(params.loc['EWMA_est'][i],ret[:t],init_step)[:,:,s:] \n", " Ht_ewma_svd [:,:,s:t] = EWMA.ewma_est_H(params.loc['EWMA_svd_est'][i],ret[:t],init_step)[:,:,s:] \n", " Ht_comp_ewma[:,:,s:t] = EWMA.ewma_est_H(params.loc['EWMA_comp_est'][i],ret[:t],init_step)[:,:,s:] \n", " Ht_comp_ewma_full[:,:,s:t] = EWMA.ewma_est_H(params.loc['EWMA_comp_est_full'][i],ret[:t],init_step)[:,:,s:]\n", " Rt_dcc[:,:,s:t] = DCC.dcceq(params.loc['DCC_est'][i],ret[:t],D_inv[:t,:],init_step)[:,:,s:] \n", " Rt_dcc_svd[:,:,s:t] = DCC.dcceq(params.loc['DCC_svd_est'][i],ret[:t],D_inv[:t,:],init_step)[:,:,s:] \n", " Rt_comp_dcc[:,:,s:t] = DCC.dcceq(params.loc['DCC_comp_est'][i],ret[:t],D_inv[:t,:],init_step)[:,:,s:] \n", " Rt_comp_dcc_full[:,:,s:t] = DCC.dcceq(params.loc['DCC_comp_est_full'][i],ret[:t],D_inv[:t,:],init_step)[:,:,s:] \n", " Ht_svech[:,:,s:t] = SV.svech_est_H(params.loc['SVECH_est'][i],ret[:t],init_step)[:,:,s:] \n", " Ht_svech_svd[:,:,s:t] = SV.svech_est_H(params.loc['SVECH_svd_est'][i],ret[:t],init_step)[:,:,s:] \n", " Ht_comp_svech[:,:,s:t]= SV.svech_est_H(params.loc['SVECH_comp_est'][i],ret[:t],init_step)[:,:,s:] \n", " Ht_comp_svech_full[:,:,s:t]= SV.svech_est_H(params.loc['SVECH_comp_est_full'][i],ret[:t],init_step)[:,:,s:] \n", " Ht_ore[:,:,s:t]= ORE.ore_est_H(params.loc['ORE_est'][i],ret[:t],init_step)[:,:,s:] \n", " Ht_ore_svd[:,:,s:t] = ORE.ore_est_H(params.loc['ORE_svd_est'][i],ret[:t],init_step)[:,:,s:] \n", " Ht_comp_ore[:,:,s:t]= ORE.ore_est_H(params.loc['ORE_comp_est'][i],ret[:t],init_step)[:,:,s:] \n", " Ht_comp_ore_full[:,:,s:t]= ORE.ore_est_H(params.loc['ORE_comp_est_full'][i],ret[:t],init_step)[:,:,s:] \n", "\n", "\n", "for t in range(T):\n", " Ht_dcc[:,:,t] = np.matmul(np.matmul(np.diagflat(D_diag[t,:]),Rt_dcc[:,:,t]),np.diagflat(D_diag[t,:]))\n", " Ht_dcc_svd[:,:,t] = np.matmul(np.matmul(np.diagflat(D_diag[t,:]),Rt_dcc_svd[:,:,t]),np.diagflat(D_diag[t,:]))\n", " Ht_comp_dcc[:,:,t] = np.matmul(np.matmul(np.diagflat(D_diag[t,:]),Rt_comp_dcc[:,:,t]),np.diagflat(D_diag[t,:]))\n", " Ht_comp_dcc_full[:,:,t] = np.matmul(np.matmul(np.diagflat(D_diag[t,:]),Rt_comp_dcc_full[:,:,t]),np.diagflat(D_diag[t,:]))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#setup matrices used in the case when the estimates with 5 day frequency for neighbouring CML were obtained\n", "Ht_comp_ewma = np.zeros((N,N,T)) \n", "Rt_comp_dcc= np.zeros((N,N,T)) \n", "Ht_comp_dcc= np.zeros((N,N,T)) \n", "Ht_comp_svech= np.zeros((N,N,T))\n", "Ht_comp_ore= np.zeros((N,N,T)) \n", "\n", "estimates_names =[\n", "'EWMA_comp_est',\n", "'DCC_comp_est',\n", "'SVECH_comp_est',\n", "'ORE_comp_est'\n", "]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "Time_cutoff =[]\n", "for t in range(300,T,5):\n", " Time_cutoff.append(t)\n", "if Time_cutoff[-1] != T:\n", " Time_cutoff.append(T)\n", "for i,t in enumerate(Time_cutoff[1:]):\n", " if t == Time_cutoff[1]:\n", " s=0\n", " else:\n", " s= Time_cutoff[i]\n", " print(i,s,t) \n", " Ht_comp_ewma[:,:,s:t] = EWMA.ewma_est_H(params.loc['EWMA',t],ret[:t],init_step)[:,:,s:] \n", " Rt_comp_dcc[:,:,s:t] = DCC.dcceq(params.loc['DCC',t],ret[:t],D_inv[:t,:],init_step)[:,:,s:] \n", " Ht_comp_svech[:,:,s:t]= SV.svech_est_H(params.loc['BEKK',t],ret[:t],init_step)[:,:,s:] \n", " Ht_comp_ore[:,:,s:t]= ORE.ore_est_H(params.loc['ORE',t],ret[:t],init_step)[:,:,s:] \n", "\n", "\n", "for t in range(T): \n", " Ht_comp_dcc[:,:,t] = np.matmul(np.matmul(np.diagflat(D_diag[t,:]),Rt_comp_dcc[:,:,t]),np.diagflat(D_diag[t,:]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculate portfolio performance\n", "a few methods of calculating and drawing portfolio performance connected variables and matrics" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "#A few methods to show porfolio performance\n", "def calc_perf(w,rets):\n", " T = w.shape[0]\n", " \n", " perf = []\n", " for t in range(T):\n", " perf_t = np.matmul(w[t,:], rets[t,:].T)\n", " perf.append(perf_t)\n", " port_perf = np.array(perf)\n", " return port_perf\n", "\n", "def calc_cum_perf(daily_ret):\n", " T = daily_ret.shape[0] \n", " daily_perf = daily_ret + 1\n", " gains= np.ones(T)\n", " \n", " for cnt in range(T):\n", " if cnt ==0 :\n", " gains[cnt] = daily_perf[cnt]\n", " else:\n", " gains[cnt] = gains[cnt-1] * daily_perf[cnt]\n", " return gains\n", "\n", "def calc_turnover(weights):\n", " T = weights.shape[0]\n", " turnover = np.zeros((T))\n", " turnover[0]=0\n", " for t in range(1,T):\n", " turnover[t] = sum(abs(weights[t,:]- weights[t-1,:]))\n", " \n", " return turnover\n", " \n", "\n", "def see_perf(ws,m_names, rets, draw_ind_stocks = True):\n", " \n", " fig, ax1 = plt.subplots(figsize=(15, 10)) \n", " # ax2 = ax1.twinx()\n", " \n", " T,N = rets.shape\n", " x= dates[300:]\n", " M= len(ws)\n", " maxgain =0 \n", " ls=\"solid\"\n", " color =[]\n", " for m,w in enumerate(ws):\n", " if m!=0:\n", " color.append(next(ax1._get_lines.prop_cycler)['color'])\n", " gains= np.ones(T)\n", " perf = calc_perf(w, rets)\n", " perf += 1\n", " \n", " for cnt,a in enumerate(perf):\n", " if cnt ==0:\n", " gains[0]= perf[0]\n", " else:\n", " gains[cnt] = gains[cnt-1] * perf[cnt]\n", " \n", " gains = np.array(gains)\n", " if maxgain < gains[-1]:\n", " maxgain = gains[-1]\n", " if draw_ind_stocks == True: \n", " ax1.plot(x,gains,lw=3, label = m_names[m])\n", " else: \n", " if ((m>5)& (m<=10)):\n", " ls= (0, (5, 1))\n", " print(m,ls,m_names[m])\n", " else:\n", " ls=\"solid\"\n", " print(m,ls,m_names[m])\n", " if m==0:\n", " col='black'\n", " else:\n", " col= color[(m-1) %5 ]\n", " ax1.plot(x,gains,lw=1.5, label = m_names[m],ls=ls, color=col)\n", " # ax2.plot(perf, color= 'blue')\n", " \n", " if draw_ind_stocks == True: \n", " for stock in range(N):\n", " ind_gain= np.ones(T)\n", " ind_gain[0] = rets[0,stock]+1\n", " for cnt in range(1,T): \n", " ind_gain[cnt] = ind_gain[cnt-1] * (rets[cnt,stock]+1)\n", " ind_gain = np.array(ind_gain)\n", " if ind_gain[-1]>maxgain:\n", " ax1.plot(ind_gain, alpha=.25, label=assets[stock] )\n", " else:\n", " ax1.plot(ind_gain, alpha=.25)\n", "\n", " ax1.legend() \n", " \n", " ax1.set_ylabel('cumulative return')\n", " # ax2.set_ylabel('daily return')\n", " \n", " return fig\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Draw port performance" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "# See properties of all models \n", "def ind_model_performance(model_names,mean_tuple,disc_mean_tuple, var_tuple,\n", " weights,shrape_tuple,turnover_tuple,\n", " eq_ret_disc=None,eq_var=None,ind =0, MA_window = 50):\n", " \n", " \n", " fig, (ax1, ax2, ax3,ax4,ax5,ax6) = plt.subplots( 6, figsize=(15, 30)) \n", " ax1.set_title('Mean portfolio delta-adjusted and non-adjusted mean returns')\n", " ax2.set_title('Delta-adjusted porfolio standard deviation')\n", " ax3.set_title('Final model weights')\n", " ax4.set_title('Sharpe ratios')\n", " ax5.set_title('Portfolio Turnover')\n", " ax6.set_title('Portfolio Turnover moving average')\n", " for i in range(len(model_names)):\n", " ax1.plot(mean_tuple[i][ind:])\n", " ax1.plot(mean_disc_tuple[i][ind:], label = model_names[i])\n", " ax2.plot(np.sqrt(var_tuple[i][ind:]), label = model_names[i])\n", " ax3.plot(weights[i,:][ind:], label = model_names[i])\n", " ax4.plot(shrape_tuple[i][ind:], label = model_names[i])\n", " ax5.plot(turnover_tuple[i][ind:], label = model_names[i])\n", " ax6.plot(pd.Series(turnover_tuple[i][ind:]).rolling(window=MA_window, min_periods=1).mean()\n", " , label = model_names[i]+' MA ' + str(MA_window))\n", "\n", " ax6.plot(turnover_tuple[i][ind:].mean()*np.ones(len(turnover_tuple[i][ind:])),lw = 2, label = model_names[i]+' mean' )\n", " \n", " \n", " if eq_ret_disc.any()!= None:\n", " ax1.plot(eq_ret_disc[ind:], label ='Stocks Equal Weights')\n", " ax2.plot(np.sqrt(eq_mean_var[ind:]), label ='Stocks Equal Weights')\n", "\n", " ax1.legend(loc= 'upper left') \n", " ax2.legend(loc= 'upper left')\n", " ax3.legend(loc= 'upper left') \n", " ax4.legend(loc= 'upper left') \n", " ax5.legend(loc= 'upper left') \n", " ax6.legend(loc= 'upper left') \n", " #fig.show\n", " \n", " return fig\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# See properties of all models \n", "def ind_model_performance2(model_names,disc_mean_tuple, var_tuple,\n", " turnover_tuple,beta_tuple,shrape_tuple,sortino_tuple,treinor_tuple,\n", " eq_tuple=None,ind =0, MA_window = 50):\n", " \n", " \n", " fig, (ax1, ax2, ax3,ax4,ax5,ax6,ax7) = plt.subplots( 7, figsize=(12, 30)) \n", " ax1.set_title('Mean portfolio returns')\n", " ax2.set_title('Portfolio standard deviation')\n", " ax3.set_title('Portfolio turnover')\n", " ax4.set_title('Porfolio beta')\n", " ax5.set_title('Sharpe ratios')\n", " ax6.set_title('Sortino ratios')\n", " ax7.set_title('Trainors index')\n", " for i in range(len(model_names)): \n", " color = next(ax1._get_lines.prop_cycler)['color']\n", " ax1.plot(mean_disc_tuple[i][ind:], label = model_names[i], color=color)\n", " ax2.plot(np.sqrt(var_tuple[i][ind:]), label = model_names[i], color=color)\n", " ax3.plot(pd.Series(turnover_tuple[i][ind:]).rolling(window=MA_window, min_periods=1).mean()\n", " ,color=color) #label = model_names[i]+' MA ' + str(MA_window)\n", " ax3.plot(turnover_tuple[i][ind:].mean()*np.ones(len(turnover_tuple[i][ind:])),lw = 2,\n", " label = model_names[i], color=color )\n", " ax4.plot(beta_tuple[i][ind:], label = model_names[i], color=color)\n", " ax5.plot(shrape_tuple[i][ind:], label = model_names[i], color=color)\n", " ax6.plot(sortino_tuple[i][ind:], label = model_names[i], color=color)\n", " ax7.plot(treinor_tuple[i][ind:], label = model_names[i], color=color)\n", "\n", "\n", " \n", " \n", " if eq_tuple!= None:\n", " lb= 'Stocks Equal Weights'\n", " color = 'black'\n", " ax1.plot(eq_tuple[0][ind:], label =lb,color=color)\n", " ax2.plot(np.sqrt(eq_tuple[1][ind:]), label =lb,color=color)\n", " \n", " ax4.plot(eq_tuple[2][ind:], label = lb,color=color)\n", " ax5.plot(eq_tuple[3][ind:], label = lb,color=color)\n", " ax6.plot(eq_tuple[4][ind:], label = lb,color=color)\n", " ax7.plot(eq_tuple[5][ind:], label = lb,color=color)\n", "\n", " ax1.legend(loc= 'upper left') \n", " ax2.legend(loc= 'upper left')\n", " ax3.legend(loc= 'upper left') \n", " ax4.legend(loc= 'upper left') \n", " ax5.legend(loc= 'upper left') \n", " ax6.legend(loc= 'upper left') \n", " ax7.legend(loc= 'upper left') \n", " #fig.show\n", " \n", " return fig" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# See properties of all models \n", "def ind_model_performance3(model_names,disc_mean_tuple, var_tuple,\n", " turnover_tuple,beta_tuple,\n", " eq_tuple=None,ind =0, MA_window = 50):\n", " \n", " linew=3 \n", "\n", " fig, (ax1, ax2, ax3,ax4) = plt.subplots( 4, figsize=(14, 22)) \n", " ax1.set_title('Mean portfolio returns')\n", " ax2.set_title('Portfolio standard devation')\n", " ax3.set_title('Portfolio turnover')\n", " ax4.set_title('Porfolio beta')\n", " # ax5.set_title('Sharpe ratios')\n", " # ax6.set_title('Sortino ratios')\n", " # ax7.set_title('Trainors index')\n", " x=dates[ind:]\n", " for i in range(len(model_names)): \n", " color = next(ax1._get_lines.prop_cycler)['color']\n", " ax1.plot(x,mean_disc_tuple[i][ind:], label = model_names[i], color=color,lw=linew)\n", " ax2.plot(x,np.sqrt(var_tuple[i][ind:]), label = model_names[i], color=color,lw=linew)\n", " ax3.plot(x,pd.Series(turnover_tuple[i][ind:]).rolling(window=MA_window, min_periods=1).mean()\n", " ,label = model_names[i],color=color,lw=linew-0.5) #label = model_names[i]+' MA ' + str(MA_window)\n", " ax3.plot(x,turnover_tuple[i][ind:].mean()*np.ones(len(turnover_tuple[i][ind:])),lw = linew+0.5,ls='dashed',\n", " color=color)\n", " ax4.plot(x,beta_tuple[i][ind:], label = model_names[i], color=color,lw=linew)\n", " # ax5.plot(shrape_tuple[i][ind:], label = model_names[i], color=color)\n", " # ax6.plot(sortino_tuple[i][ind:], label = model_names[i], color=color)\n", " # ax7.plot(treinor_tuple[i][ind:], label = model_names[i], color=color)\n", "\n", "\n", " \n", " \n", " if eq_tuple!= None:\n", " lb= 'Stocks Equal Weights'\n", " color = 'black'\n", " ax1.plot(x,eq_tuple[0][ind:], label =lb,color=color,lw=linew)\n", " ax2.plot(x,np.sqrt(eq_tuple[1][ind:]), label =lb,color=color,lw=linew)\n", " \n", " ax4.plot(x,eq_tuple[2][ind:], label = lb,color=color,lw=linew)\n", " # ax5.plot(eq_tuple[3][ind:], label = lb,color=color)\n", " # ax6.plot(eq_tuple[4][ind:], label = lb,color=color)\n", " # ax7.plot(eq_tuple[5][ind:], label = lb,color=color)\n", "\n", " ax1.legend(loc= 'upper left') \n", " ax2.legend(loc= 'upper left')\n", " ax3.legend(loc= 'upper left') \n", " ax4.legend(loc= 'upper left') \n", " # ax5.legend(loc= 'upper left') \n", " # ax6.legend(loc= 'upper left') \n", " # ax7.legend(loc= 'upper left') \n", " #fig.show\n", " \n", " return fig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mean-Variance problem under forecast combination" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "#Two methods to estimate stock expected returns\n", "def rolling_mean(mat, wind): #get column-wise rolling mean \n", " rmean= []\n", " T,N= ret.shape\n", " for col in range(N):\n", " rmean.append( pd.Series(mat[:,col]).rolling(window=wind, min_periods=1).mean())\n", " \n", " mu = np.array(rmean).T\n", " return mu \n", "\n", "def discounted_rolling_mean(mat,deltas):\n", " dmean= []\n", " T,N= ret.shape\n", " d = deltas[1]\n", " for t in range(T):\n", " if d == 1:\n", " dmean.append(1/(t+1)*np.matmul(deltas[:t+1][::-1].T, mat[:t+1]))\n", " else:\n", " dmean.append((1 - d)/(1 - d**(max(t,2)-1))*np.matmul(deltas[:t+1][::-1].T, mat[:t+1]))\n", " \n", " return np.array(dmean)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "def meanvar(H, mu,mu_0):\n", " w =[]\n", " T = H.shape[2]\n", " for t in range(T):\n", " Ht=H[:,:,t]\n", " if t == 0:\n", " mu_t = np.array(mu[t,:],ndmin = 2).T\n", " else:\n", " mu_t = np.array(mu[t-1,:],ndmin = 2).T\n", " \n", " Ht_inv = np.linalg.inv(Ht) # u had **2 here for some reason\n", " ones = np.array(np.ones(N), ndmin=2).T #column vector of ones\n", " A = np.matmul(np.matmul(mu_t.T, Ht_inv), mu_t)\n", " B = np.matmul(np.matmul(mu_t.T, Ht_inv), ones)\n", " C = np.matmul(np.matmul(ones.T,Ht_inv),ones)\n", " \n", " nom = np.matmul(mu_t, mu_0*C -B) + np.matmul(ones , A - mu_0 * B) \n", " denom = A * C - B**2\n", "\n", " w.append(np.divide(np.matmul(Ht_inv, nom), denom))\n", " \n", " weights = np.array(w)\n", " weights = np.squeeze(weights)\n", " return weights" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Min var problem under forecast information" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "def minvar(H):\n", " w =[]\n", " T = H.shape[2]\n", " for t in range(T):\n", " Ht=H[:,:,t]\n", " Ht_inv = np.linalg.inv(Ht)\n", " ones = np.array(np.ones(N), ndmin=2).T #column vector of ones\n", " w_t = np.matmul(ones.T,Ht_inv) #not normalized\n", " denom = np.matmul(w_t, ones)\n", " w_t = np.divide(w_t.T, denom)\n", " w.append(w_t)\n", " \n", " weights = np.array(w)\n", " weights = np.squeeze(weights)\n", " return weights" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Risk adjusted performance metrics" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "def port_return(w, ret, deltas):\n", " d = deltas[1]\n", " T = ret.shape[0]\n", " day_perf =np.zeros(T)\n", " day_perf_disc= np.zeros(T)\n", " mean_port_ret= np.zeros(T)\n", " mean_port_ret_disc= np.zeros(T)\n", " \n", " for t in range(0,T):\n", " day_perf[t] = np.array(np.matmul(np.array(ret[t,:],ndmin =2), w[t,:] ),ndmin =2)\n", " mean_port_ret[t] = 1/(max(t,2)-1)*sum(day_perf[:t])\n", " if d!=1:\n", " mean_port_ret_disc[t]= (1 - d)/(1 - d**(max(t,2)-1))*np.matmul(day_perf[:t+1],deltas[:t+1][::-1])\n", " else:\n", " mean_port_ret_disc[t]= 1/max(t+1,1)*np.matmul(day_perf[:t+1],deltas[:t+1][::-1])\n", " #[::-1] reverses the array\n", " \n", " return day_perf, mean_port_ret, mean_port_ret_disc\n", "\n", "def port_return_costadj(w, ret,turn):\n", " \n", " T = ret.shape[0]\n", " day_perf_costadj =np.zeros(T)\n", " c_list =np.array([0,5,10,15,20])*0.0001\n", " gain=[1,1,1,1,1]\n", " \n", " for t in range(0,T):\n", " day_perf_costadj[t] = np.array(np.matmul(np.array(ret[t,:]+1,ndmin =2), w[t,:] ),ndmin =2)\n", " for i,c in enumerate(c_list):\n", " gain[i] = gain[i]*day_perf_costadj[t]*(1-c*turn[t])\n", " return gain\n", "\n", "\n", "def port_variance(daily_perf, mu_m, deltas,start=0):\n", " T = daily_perf.shape[0]\n", " d= deltas[1]\n", " m_var= np.zeros(T)\n", " for t in range(start,T):\n", " if d!=1:\n", " m_var[t] = np.array(np.matmul(((daily_perf[:t+1]-mu_m[t])**2), deltas[:t+1][::-1]).sum()*((1 - d)/(1 - d**(max(t,2)-1))))\n", " else:\n", " m_var[t] = np.array(np.matmul(((daily_perf[:t+1]-mu_m[t])**2), deltas[:t+1][::-1]).sum()*1/t)\n", " return m_var\n", "\n", "def port_downside_variance(daily_perf, mu_m, deltas,start=0):\n", " T = daily_perf.shape[0]\n", " d= deltas[1]\n", " m_down_var= np.zeros(T)\n", " down_perf= np.where(daily_perf>0, 0, daily_perf)\n", " for t in range(start,T):\n", " if d!=1:\n", " m_down_var[t] = np.array(np.matmul((down_perf[:t+1]**2), deltas[:t+1][::-1]).sum()*((1 - d)/(1 - d**(max(t,2)-1))))\n", " else:\n", " m_down_var[t] = np.array(np.matmul((down_perf[:t+1]**2), deltas[:t+1][::-1]).sum()*1/t)\n", " \n", " m_down_var= np.where(((m_down_var<0.000001)&(m_down_var>0)), 0.000001, m_down_var)\n", " m_down_var= np.where(((m_down_var>-0.000001)&(m_down_var<0)), -0.000001, m_down_var)\n", " return m_down_var\n", " \n", "def sharpe_ratio(model_ret,model_variance,minvar=False):\n", " if minvar==True:\n", " return np.ones(model_ret.shape[0])/np.sqrt(model_variance)\n", " \n", " return (model_ret/np.sqrt(model_variance))\n", "\n", "#TODO\n", "def sortino_ratio(model_ret,model_down_variance,minvar=False):\n", " if minvar == True:\n", " return np.ones(model_ret.shape[0])/np.sqrt(model_down_variance)\n", " SR = model_ret/np.sqrt(model_down_variance)\n", " return SR\n", " \n", "def port_beta(daily_perf,mu_m, index_dperf,deltas,start =1):\n", " T = daily_perf.shape[0]\n", " d= deltas[1]\n", " betas= np.zeros(T)\n", " mean_index_ret = np.zeros(T)\n", " betas[0]=1\n", " for t in range(start,T):\n", " mean_index_ret[t] = 1/(max(t,2)-1)*sum(index_dperf[:min(t+1,T)])\n", " if d!=1:\n", " betas[t] = np.array(np.matmul(((daily_perf[:t+1] - mu_m[t])*(index_dperf[:t+1,0] - mean_index_ret[t])).T, deltas[:t+1][::-1]).sum()*((1 - d)/(1 - d**(max(t,2)-1))))\n", " denom = np.matmul(((index_dperf[:t+1,0] - mean_index_ret[t])**2).T, deltas[:t+1][::-1]).sum()*((1 - d)/(1 - d**(max(t,2)-1)))\n", " else:\n", " betas[t] = np.array(np.matmul(((daily_perf[:t+1] - mu_m[t])*(index_dperf[:t+1,0] - mean_index_ret[t])).T, deltas[:t+1][::-1]).sum()*1/t)\n", " denom = np.matmul(((index_dperf[:t+1,0] - mean_index_ret[t])**2).T, deltas[:t+1][::-1]).sum()*1/t\n", " betas[t]= betas[t]/denom\n", " betas= np.where(((betas<0.01)&(betas>0)), 0.01, betas)\n", " betas= np.where(((betas>-0.01)&(betas<0)), -0.01, betas)\n", " return betas\n", "\n", "def trainors_index(model_ret,beta): \n", " TI = model_ret/ abs(beta)\n", " return TI\n", " \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduce a discount parameter" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for delta in (0.99,0.95,0.9,0.8,0.7):\n", " deltas = np.zeros(T)\n", " deltas[0] = 1\n", " deltas[1] = delta \n", " for t in range(2,T):\n", " deltas[t] = deltas[t-1]*delta\n", " plt.plot(deltas[:100], label = r\"$\\delta=$\" + str(delta))\n", " plt.legend() \n", " plt.xlabel(r'$t $' , color=\"black\", fontsize = 15)\n", " plt.ylabel(r'$\\delta^t $' , color=\"black\", fontsize = 15)\n", " #plt.savefig('dicountParam.pdf')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "delta= 1\n", "deltas = np.zeros(T)\n", "deltas[0] = 1\n", "deltas[1] = delta\n", "for t in range(2,T):\n", " deltas[t] = deltas[t-1]*delta\n", "plt.plot(deltas[:200])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Estimating the combined matix\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mean-var " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "use_composite_est = True" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mu = rolling_mean(ret, 500) #column vector of expected returns\n", "#mu = ret.mean(axis = 0)\n", "#mu = discounted_rolling_mean(ret,deltas) #Smaller turnover when delta relatively close to 1\n", "\n", "daily_exp_ret = 0.0005\n", "\n", "if use_composite_est==False:\n", " mean_w_ccc = meanvar(Ht_ccc, mu, daily_exp_ret)\n", " mean_w_svech = meanvar(Ht_svech, mu, daily_exp_ret)\n", " mean_w_ore = meanvar(Ht_ore, mu, daily_exp_ret)\n", " mean_w_dcc = meanvar(Ht_dcc, mu, daily_exp_ret)\n", " mean_w_ewma = meanvar(Ht_ewma, mu, daily_exp_ret)\n", "else:\n", " mean_w_svech = meanvar(Ht_comp_svech, mu, daily_exp_ret)\n", " mean_w_ore = meanvar(Ht_comp_ore, mu, daily_exp_ret)\n", " mean_w_dcc = meanvar(Ht_comp_dcc, mu, daily_exp_ret)\n", " mean_w_ewma = meanvar(Ht_comp_ewma, mu, daily_exp_ret)\n", " mean_w_ccc = meanvar(Ht_ccc, mu, daily_exp_ret)\n", " \n", "equal_weights = np.ones((T,N))*1/N\n", "\n", "mean_turn_svech = calc_turnover(mean_w_svech)\n", "mean_turn_ore = calc_turnover(mean_w_ore)\n", "mean_turn_dcc = calc_turnover(mean_w_dcc)\n", "mean_turn_ewma = calc_turnover(mean_w_ewma)\n", "mean_turn_ccc = calc_turnover(mean_w_ccc)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "svech_mean_dperf, svech_mean_ret, svech_mean_ret_disc = port_return(mean_w_svech,ret,deltas) \n", "ore_mean_dperf, ore_mean_ret, ore_mean_ret_disc = port_return(mean_w_ore,ret,deltas)\n", "ccc_mean_dperf, ccc_mean_ret, ccc_mean_ret_disc = port_return(mean_w_ccc,ret,deltas)\n", "dcc_mean_dperf, dcc_mean_ret, dcc_mean_ret_disc = port_return(mean_w_dcc,ret,deltas)\n", "ewma_mean_dperf, ewma_mean_ret, ewma_mean_ret_disc = port_return(mean_w_ewma,ret,deltas)\n", "eq_mean_dperf, eq_mean_ret, eq_mean_ret_disc = port_return(equal_weights,ret,deltas)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#use undiscounted mean performance here? \n", "svech_mean_var = port_variance(svech_mean_dperf,svech_mean_ret,deltas)\n", "ore_mean_var = port_variance(ore_mean_dperf,ore_mean_ret,deltas)\n", "dcc_mean_var = port_variance(dcc_mean_dperf,dcc_mean_ret,deltas)\n", "ccc_mean_var = port_variance(ccc_mean_dperf,ccc_mean_ret,deltas)\n", "ewma_mean_var = port_variance(ewma_mean_dperf,ewma_mean_ret,deltas)\n", "eq_mean_var = port_variance(eq_mean_dperf,eq_mean_ret,deltas)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "svech_mean_downvar = port_downside_variance(svech_mean_dperf,svech_mean_ret,deltas)\n", "ore_mean_downvar = port_downside_variance(ore_mean_dperf,ore_mean_ret,deltas)\n", "dcc_mean_downvar = port_downside_variance(dcc_mean_dperf,dcc_mean_ret,deltas)\n", "ccc_mean_downvar = port_downside_variance(ccc_mean_dperf,ccc_mean_ret,deltas)\n", "ewma_mean_downvar = port_downside_variance(ewma_mean_dperf,ewma_mean_ret,deltas)\n", "eq_mean_downvar = port_downside_variance(eq_mean_dperf,eq_mean_ret,deltas)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "svech_mean_beta = port_beta(svech_mean_dperf, svech_mean_ret,SP_ret,deltas)\n", "ore_mean_beta = port_beta(ore_mean_dperf, ore_mean_ret,SP_ret,deltas)\n", "ccc_mean_beta = port_beta(ccc_mean_dperf, ccc_mean_ret,SP_ret,deltas)\n", "dcc_mean_beta = port_beta(dcc_mean_dperf, dcc_mean_ret,SP_ret,deltas)\n", "ewma_mean_beta = port_beta(ewma_mean_dperf, ewma_mean_ret,SP_ret,deltas) \n", "eq_mean_beta = port_beta(eq_mean_dperf, eq_mean_ret,SP_ret,deltas) " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "svech_mean_SR = sharpe_ratio(svech_mean_ret_disc, svech_mean_var)\n", "ore_mean_SR = sharpe_ratio(ore_mean_ret_disc,ore_mean_var)\n", "ccc_mean_SR = sharpe_ratio(ccc_mean_ret_disc,ccc_mean_var)\n", "dcc_mean_SR = sharpe_ratio(dcc_mean_ret_disc,dcc_mean_var)\n", "ewma_mean_SR = sharpe_ratio(ewma_mean_ret_disc, ewma_mean_var)\n", "eq_mean_SR = sharpe_ratio(eq_mean_ret_disc, eq_mean_var)\n", "\n", "svech_mean_SoR = sortino_ratio(svech_mean_ret_disc, svech_mean_downvar)\n", "ore_mean_SoR = sortino_ratio(ore_mean_ret_disc,ore_mean_downvar)\n", "ccc_mean_SoR = sortino_ratio(ccc_mean_ret_disc,ccc_mean_downvar)\n", "dcc_mean_SoR = sortino_ratio(dcc_mean_ret_disc,dcc_mean_downvar)\n", "ewma_mean_SoR = sortino_ratio(ewma_mean_ret_disc, ewma_mean_downvar)\n", "eq_mean_SoR = sortino_ratio(eq_mean_ret_disc, eq_mean_downvar)\n", "\n", "\n", "svech_mean_TI = trainors_index(svech_mean_ret_disc, svech_mean_beta)\n", "ore_mean_TI = trainors_index(ore_mean_ret_disc, ore_mean_beta)\n", "ccc_mean_TI = trainors_index(ccc_mean_ret_disc, ccc_mean_beta)\n", "dcc_mean_TI = trainors_index(dcc_mean_ret_disc, dcc_mean_beta)\n", "ewma_mean_TI = trainors_index(ewma_mean_ret_disc, ewma_mean_beta)\n", "eq_mean_TI = trainors_index(eq_mean_ret_disc, eq_mean_beta)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#keep im mind that weights at time t have to detrmine positions for the next day t+1\n", "ni = 1 #model weights aplification\n", "\n", "models = (svech_mean_SR,ore_mean_SR,ccc_mean_SR,dcc_mean_SR,ewma_mean_SR)\n", "model_names =('BEKK','ORE','CCC','DCC','EWMA')\n", "M = len(models)\n", "init_positions = np.ones(M) * 1/M\n", "\n", "stack = np.stack(arrays = models)\n", "stack= np.where(stack<0, 0, stack)\n", "stack = stack **ni\n", "\n", "cum_weights = stack.sum(axis = 0 )\n", "model_mean_weights = stack / cum_weights\n", "model_mean_weights[:,0] = init_positions\n", "\n", "indx= np.where(cum_weights<=0) #cases when all SRs are below zero - choose equal weights\n", "model_mean_weights[:,indx] = 1/M\n", "model_mean_weights[:,0] = 1/len(models) # initial setup\n", "\n", "shifted_mean_weights = np.roll(model_mean_weights, -1, axis=1)\n", "\n", "shifted_mean_weights[:,0] = init_positions\n", "\n", "shifted_mean_weights = np.insert(model_mean_weights, 0, init_positions, axis=1) # initial setup\n", "shifted_mean_weights = np.delete(shifted_mean_weights,obj = -1, axis =1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "ind =300\n", "mean_disc_tuple = (svech_mean_ret_disc,ore_mean_ret_disc,ccc_mean_ret_disc,\n", " dcc_mean_ret_disc,ewma_mean_ret_disc,eq_mean_ret_disc)\n", "mean_tuple = (svech_mean_ret,ore_mean_ret,ccc_mean_ret,dcc_mean_ret,\n", " ewma_mean_ret,eq_mean_ret)\n", "var_mean_tuple = (svech_mean_var, ore_mean_var,ccc_mean_var, dcc_mean_var,ewma_mean_var, eq_mean_var)\n", "mean_sharpe_tuple = (svech_mean_SR,ore_mean_SR,ccc_mean_SR,dcc_mean_SR,ewma_mean_SR)\n", "mean_turnover_tuple = (mean_turn_svech, mean_turn_ore,mean_turn_ccc, mean_turn_dcc, mean_turn_ewma)\n", "\n", "ind_model_performance(model_names,mean_tuple,mean_disc_tuple,var_mean_tuple,\n", " shifted_mean_weights,mean_sharpe_tuple,mean_turnover_tuple,\n", " eq_mean_ret_disc,eq_mean_var,init_step,25)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "ind =init_step\n", "mean_disc_tuple = (svech_mean_ret_disc,ore_mean_ret_disc,ccc_mean_ret_disc,\n", " dcc_mean_ret_disc,ewma_mean_ret_disc,eq_mean_ret_disc)\n", "mean_tuple = (svech_mean_ret,ore_mean_ret,ccc_mean_ret,dcc_mean_ret,\n", " ewma_mean_ret,eq_mean_ret)\n", "var_mean_tuple = (svech_mean_var, ore_mean_var,ccc_mean_var, dcc_mean_var,ewma_mean_var, eq_mean_var)\n", "mean_turnover_tuple = (mean_turn_svech, mean_turn_ore,mean_turn_ccc, mean_turn_dcc, mean_turn_ewma)\n", "mean_beta_tuple= (svech_mean_beta,ore_mean_beta,ccc_mean_beta,dcc_mean_beta,ewma_mean_beta)\n", "\n", "mean_sharpe_tuple = (svech_mean_SR,ore_mean_SR,ccc_mean_SR,dcc_mean_SR,ewma_mean_SR)\n", "mean_sortino_tuple = (svech_mean_SoR,ore_mean_SoR,ccc_mean_SoR,dcc_mean_SoR,ewma_mean_SoR)\n", "mean_trainor_tuple = (svech_mean_TI,ore_mean_TI,ccc_mean_TI,dcc_mean_TI,ewma_mean_TI)\n", "\n", "eqw_tuple= (eq_mean_ret_disc,eq_mean_var,eq_mean_beta,eq_mean_SR,eq_mean_SoR,eq_mean_TI)\n", "\n", "figure = ind_model_performance2(model_names,mean_disc_tuple,var_mean_tuple,\n", " mean_turnover_tuple,mean_beta_tuple ,mean_sharpe_tuple,mean_sharpe_tuple,\n", " mean_trainor_tuple,eqw_tuple,300,25)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#figure.savefig(\"Meanvar_candm1.pdf\")\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ind =init_step\n", "mean_disc_tuple = (svech_mean_ret_disc,ore_mean_ret_disc,ccc_mean_ret_disc,\n", " dcc_mean_ret_disc,ewma_mean_ret_disc,eq_mean_ret_disc)\n", "mean_tuple = (svech_mean_ret,ore_mean_ret,ccc_mean_ret,dcc_mean_ret,\n", " ewma_mean_ret,eq_mean_ret)\n", "var_mean_tuple = (svech_mean_var, ore_mean_var,ccc_mean_var, dcc_mean_var,ewma_mean_var, eq_mean_var)\n", "mean_turnover_tuple = (mean_turn_svech, mean_turn_ore,mean_turn_ccc, mean_turn_dcc, mean_turn_ewma)\n", "mean_beta_tuple= (svech_mean_beta,ore_mean_beta,ccc_mean_beta,dcc_mean_beta,ewma_mean_beta)\n", "\n", "mean_sharpe_tuple = (svech_mean_SR,ore_mean_SR,ccc_mean_SR,dcc_mean_SR,ewma_mean_SR)\n", "mean_sortino_tuple = (svech_mean_SoR,ore_mean_SoR,ccc_mean_SoR,dcc_mean_SoR,ewma_mean_SoR)\n", "mean_trainor_tuple = (svech_mean_TI,ore_mean_TI,ccc_mean_TI,dcc_mean_TI,ewma_mean_TI)\n", "\n", "eqw_tuple= (eq_mean_ret_disc,eq_mean_var,eq_mean_beta)\n", "\n", "figure = ind_model_performance3(model_names,mean_disc_tuple,var_mean_tuple,\n", " mean_turnover_tuple,mean_beta_tuple ,eqw_tuple,300,25)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "H_mean_comb = np.zeros((N,N,T))\n", "H_mean_eqw_comb = np.zeros((N,N,T))\n", "M=shifted_mean_weights.shape[0]\n", "\n", "if use_composite_est==False:\n", " for t in range(T):\n", " H_mean_comb[:,:,t] = (shifted_mean_weights[0,t]* Ht_svech[:,:,t] + shifted_mean_weights[1,t]* Ht_ore[:,:,t] + \\\n", " shifted_mean_weights[2,t] * Ht_ccc[:,:,t]+shifted_mean_weights[3,t] * Ht_dcc[:,:,t] + \\\n", " shifted_mean_weights[4,t]* Ht_ewma[:,:,t])\n", "\n", " H_mean_eqw_comb[:,:,t] = (1/M* Ht_svech[:,:,t] + 1/M* Ht_ore[:,:,t] + \n", " 1/M * Ht_ccc[:,:,t] +1/M * Ht_dcc[:,:,t] + 1/M* Ht_ewma[:,:,t])\n", "else:\n", " for t in range(T):\n", " H_mean_comb[:,:,t] = (shifted_mean_weights[0,t]* Ht_comp_svech[:,:,t] + shifted_mean_weights[1,t]* Ht_comp_ore[:,:,t]+ \\\n", " shifted_mean_weights[2,t] * Ht_ccc[:,:,t]+shifted_mean_weights[3,t] * Ht_comp_dcc[:,:,t] + \\\n", " shifted_mean_weights[4,t]* Ht_comp_ewma[:,:,t])\n", "\n", " H_mean_eqw_comb[:,:,t] = (1/M* Ht_comp_svech[:,:,t] + 1/M* Ht_comp_ore[:,:,t] + \n", " 1/M * Ht_ccc[:,:,t] +1/M * Ht_comp_dcc[:,:,t] + 1/M* Ht_comp_ewma[:,:,t])\n", " \n", "\n", "mean_final_w = meanvar(H_mean_comb, mu, daily_exp_ret)\n", "mean_final_eq_w = meanvar(H_mean_eqw_comb, mu, daily_exp_ret)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "a = 1+calc_perf(mean_final_w,ret)\n", "gain=1\n", "for i in range(len(a)):\n", " gain = a[i]*gain\n", "gain" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "cut_start = valid_wind\n", "cut_end = T\n", "M_consid = 10\n", "m_weights= (mean_final_w[cut_start:cut_end],mean_w_svech[cut_start:cut_end], mean_w_ore[cut_start:cut_end],\n", " mean_w_ccc[cut_start:cut_end],mean_w_dcc[cut_start:cut_end], mean_w_ewma[cut_start:cut_end],\n", " equal_weights[cut_start:cut_end],mean_final_eq_w[cut_start:cut_end] )\n", "m_names = ('COMBINATION','SVECH', 'ORE','CCC','DCC','EWMA','EQUAL WEIGHTS STOCKS','EQUAL WEIGHTS MODELS')\n", "\n", "\n", "see_perf(m_weights[:M_consid],m_names[:M_consid],ret[cut_start:cut_end], False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Min var" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "use_composite_est = True" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# model uses information from t-1 to get the estimate of H_t\n", "# the estimate H_t is used to obtain weights to be used at trading day t\n", "if use_composite_est==False:\n", " min_w_svech = minvar(Ht_svech)\n", " min_w_ore = minvar(Ht_ore)\n", " min_w_ccc = minvar(Ht_ccc)\n", " min_w_dcc = minvar(Ht_dcc)\n", " min_w_ewma = minvar(Ht_ewma)\n", "else:\n", " min_w_svech = minvar(Ht_comp_svech)\n", " min_w_ore = minvar(Ht_comp_ore)\n", " min_w_ccc = minvar(Ht_ccc)\n", " min_w_dcc = minvar(Ht_comp_dcc)\n", " min_w_ewma = minvar(Ht_comp_ewma)\n", "equal_weights = np.ones((T,N))*1/N\n", "\n", "\n", "min_turn_svech = calc_turnover(min_w_svech)\n", "min_turn_ore = calc_turnover(min_w_ore)\n", "min_turn_ccc = calc_turnover(min_w_ccc)\n", "min_turn_dcc = calc_turnover(min_w_dcc)\n", "min_turn_ewma = calc_turnover(min_w_ewma)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#The weights for day t estimated using info available at t-1 are used\n", "#To get the real performance during the trading day t using stocks returns observed at time t\n", "svech_min_dperf, svech_min_ret, svech_min_ret_disc = port_return(min_w_svech,ret,deltas) \n", "ore_min_dperf, ore_min_ret, ore_min_ret_disc = port_return(min_w_ore,ret,deltas)\n", "ccc_min_dperf, ccc_min_ret, ccc_min_ret_disc = port_return(min_w_ccc,ret,deltas)\n", "dcc_min_dperf, dcc_min_ret, dcc_min_ret_disc = port_return(min_w_dcc,ret,deltas)\n", "ewma_min_dperf, ewma_min_ret, ewma_min_ret_disc = port_return(min_w_ewma,ret,deltas)\n", "eq_min_dperf, eq_min_ret, eq_min_ret_disc = port_return(equal_weights,ret,deltas)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#use undiscounted mean performance here? Perhpas use yeally rolling mean here instead\n", "#info from time t is used to estimate profolio performance during the trading day t\n", "svech_min_var = port_variance(svech_min_dperf,svech_min_ret,deltas)\n", "ore_min_var = port_variance(ore_min_dperf,ore_min_ret,deltas)\n", "ccc_min_var = port_variance(ccc_min_dperf,ccc_min_ret,deltas)\n", "dcc_min_var = port_variance(dcc_min_dperf,dcc_min_ret,deltas)\n", "ewma_min_var = port_variance(ewma_min_dperf,ewma_min_ret,deltas)\n", "eq_min_var = port_variance(eq_min_dperf,eq_min_ret,deltas)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "svech_min_downvar = port_downside_variance(svech_min_dperf,svech_min_ret,deltas)\n", "ore_min_downvar = port_downside_variance(ore_min_dperf,ore_min_ret,deltas)\n", "dcc_min_downvar = port_downside_variance(dcc_min_dperf,dcc_min_ret,deltas)\n", "ccc_min_downvar = port_downside_variance(ccc_min_dperf,ccc_min_ret,deltas)\n", "ewma_min_downvar = port_downside_variance(ewma_min_dperf,ewma_min_ret,deltas)\n", "eq_min_downvar = port_downside_variance(eq_min_dperf,eq_min_ret,deltas)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "svech_min_beta = port_beta(svech_min_dperf, svech_min_ret,SP_ret,deltas)\n", "ore_min_beta = port_beta(ore_min_dperf, ore_min_ret,SP_ret,deltas)\n", "ccc_min_beta = port_beta(ccc_min_dperf, ccc_min_ret,SP_ret,deltas)\n", "dcc_min_beta = port_beta(dcc_min_dperf, dcc_min_ret,SP_ret,deltas)\n", "ewma_min_beta = port_beta(ewma_min_dperf, ewma_min_ret,SP_ret,deltas) \n", "eq_min_beta = port_beta(eq_min_dperf, eq_min_ret,SP_ret,deltas)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#Sharpe ratios of the trading day t are obtained from information available at the end of trading day t\n", "svech_min_SR = sharpe_ratio(svech_min_ret_disc, svech_min_var)\n", "ore_min_SR = sharpe_ratio(ore_min_ret_disc,ore_min_var)\n", "ccc_min_SR = sharpe_ratio(ccc_min_ret_disc,ccc_min_var)\n", "dcc_min_SR = sharpe_ratio(dcc_min_ret_disc,dcc_min_var)\n", "ewma_min_SR = sharpe_ratio(ewma_min_ret_disc, ewma_min_var)\n", "eq_min_SR = sharpe_ratio(eq_min_ret_disc, eq_min_var)\n", "\n", "svech_min_SoR = sortino_ratio(svech_min_ret_disc, svech_min_downvar)\n", "ore_min_SoR = sortino_ratio(ore_min_ret_disc,ore_min_downvar)\n", "ccc_min_SoR = sortino_ratio(ccc_min_ret_disc,ccc_min_downvar)\n", "dcc_min_SoR = sortino_ratio(dcc_min_ret_disc,dcc_min_downvar)\n", "ewma_min_SoR = sortino_ratio(ewma_min_ret_disc, ewma_min_downvar)\n", "eq_min_SoR = sortino_ratio(eq_min_ret_disc, eq_min_downvar)\n", "\n", "\n", "svech_min_TI = trainors_index(svech_min_ret_disc, svech_min_beta)\n", "ore_min_TI = trainors_index(ore_min_ret_disc, ore_min_beta)\n", "ccc_min_TI = trainors_index(ccc_min_ret_disc, ccc_min_beta)\n", "dcc_min_TI = trainors_index(dcc_min_ret_disc, dcc_min_beta)\n", "ewma_min_TI = trainors_index(ewma_min_ret_disc, ewma_min_beta)\n", "eq_min_TI = trainors_index(eq_min_ret_disc, eq_min_beta)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#keep im mind that weights at time t have to detrmine positions for the next day t+1\n", "#To get the weights for trading day t, we need to use the sharpe ratios available at the end of trading day t-1\n", "ni = 1\n", "models = (svech_min_SR,ore_min_SR,ccc_min_SR,dcc_min_SR,ewma_min_SR)\n", "model_names = ('svech_min','ore_min','ccc_min','dcc_min','ewma_min')\n", "M = len(models)\n", "init_positions = np.ones(M) * 1/M\n", "\n", "stack = np.stack(arrays = models)\n", "stack= np.where(stack<0, 0, stack)\n", "stack = stack**ni\n", "cum_weights = stack.sum(axis = 0 )\n", "model_min_weights = stack / cum_weights\n", "model_min_weights[:,0] = init_positions\n", "\n", "\n", "indx= np.where(cum_weights<=0) #cases when all SRs are below zero - choose equal weights\n", "model_min_weights[:,indx] = 1/M\n", "\n", "shifted_min_weights = np.roll(model_min_weights, -1, axis=1)\n", "\n", "#make weights at time t dependent only on info from t-1 by shifting them\n", "shifted_min_weights[:,0] = init_positions\n", "print(shifted_min_weights.shape)\n", "shifted_min_weights = np.insert(model_min_weights, 0, init_positions, axis=1) # initial setup\n", "shifted_min_weights = np.delete(shifted_min_weights,obj = -1, axis =1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "min_avg_tuple = (svech_min_ret,ore_min_ret,ccc_min_ret,dcc_min_ret,ewma_min_ret,\n", " eq_min_ret)\n", "min_disc_mean_tuple = (svech_min_ret_disc,ore_min_ret_disc,ccc_min_ret_disc,dcc_min_ret_disc,\n", " ewma_min_ret_disc,eq_min_ret_disc)\n", "min_var_tuple = (svech_min_var, ore_min_var,ccc_min_var, dcc_min_var,ewma_min_var,eq_min_var)\n", "min_sharpe_tuple = (svech_min_SR, ore_min_SR,ccc_min_SR, dcc_min_SR,ewma_min_SR,eq_min_SR)\n", "min_sortino_tuple = (svech_min_SoR, ore_min_SoR,ccc_min_SoR, dcc_min_SoR,ewma_min_SoR,eq_min_SoR)\n", "min_sharpe_tuple = (svech_min_SR, ore_min_SR,ccc_min_SR, dcc_min_SR,ewma_min_SR,eq_min_SR)\n", "min_turnover_tuple = (min_turn_svech, min_turn_ore,min_turn_ccc,min_turn_dcc, min_turn_ewma)\n", "\n", "ind_model_performance(model_names,min_avg_tuple,min_disc_mean_tuple,min_var_tuple,\n", " shifted_min_weights,min_sharpe_tuple,min_turnover_tuple,\n", " eq_min_ret_disc,eq_min_var,init_step,25)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "ind =init_step\n", "mean_disc_tuple = (svech_min_ret_disc,ore_min_ret_disc,ccc_min_ret_disc,\n", " dcc_min_ret_disc,ewma_min_ret_disc,eq_min_ret_disc)\n", "mean_tuple = (svech_min_ret,ore_min_ret,ccc_min_ret,dcc_min_ret,\n", " ewma_min_ret,eq_min_ret)\n", "var_min_tuple = (svech_min_var, ore_min_var,ccc_min_var, dcc_min_var,ewma_min_var, eq_min_var)\n", "mean_turnover_tuple = (mean_turn_svech, mean_turn_ore,mean_turn_ccc, mean_turn_dcc, mean_turn_ewma)\n", "mean_beta_tuple= (svech_min_beta,ore_min_beta,ccc_min_beta,dcc_min_beta,ewma_min_beta)\n", "\n", "mean_sharpe_tuple = (svech_min_SR,ore_min_SR,ccc_min_SR,dcc_min_SR,ewma_min_SR)\n", "mean_sortino_tuple = (svech_min_SoR,ore_min_SoR,ccc_min_SoR,dcc_min_SoR,ewma_min_SoR)\n", "mean_trainor_tuple = (svech_min_TI,ore_min_TI,ccc_min_TI,dcc_min_TI,ewma_min_TI)\n", "\n", "eqw_tuple= (eq_min_ret_disc,eq_min_var,eq_min_beta)\n", "\n", "figure = ind_model_performance3(model_names,mean_disc_tuple,var_min_tuple,\n", " mean_turnover_tuple,mean_beta_tuple ,eqw_tuple,300,25)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "H_min_comb = np.zeros((N,N,T))\n", "H_min_eqw_comb = np.zeros((N,N,T))\n", "M=shifted_mean_weights.shape[0]\n", "\n", "if use_composite_est==False:\n", " for t in range(T):\n", " H_min_comb[:,:,t] = (shifted_min_weights[0,t]* Ht_svech[:,:,t] + shifted_min_weights[1,t]* Ht_ore[:,:,t] + \n", " shifted_min_weights[2,t] * Ht_ccc[:,:,t]+ shifted_min_weights[3,t] * Ht_dcc[:,:,t] + shifted_min_weights[4,t]* Ht_ewma[:,:,t])\n", "\n", " H_min_eqw_comb[:,:,t] = (1/M* Ht_svech[:,:,t] + 1/M* Ht_ore[:,:,t] + \n", " 1/M * Ht_ccc[:,:,t] + 1/M * Ht_dcc[:,:,t] + 1/M* Ht_ewma[:,:,t])\n", "else:\n", " for t in range(T):\n", " H_min_comb[:,:,t] = (shifted_min_weights[0,t]* Ht_comp_svech[:,:,t] + shifted_min_weights[1,t]* Ht_comp_ore[:,:,t] + \n", " shifted_min_weights[2,t] * Ht_ccc[:,:,t]+ shifted_min_weights[3,t] * Ht_comp_dcc[:,:,t] + shifted_min_weights[4,t]* Ht_comp_ewma[:,:,t])\n", "\n", " H_min_eqw_comb[:,:,t] = (1/M* Ht_comp_svech[:,:,t] + 1/M* Ht_comp_ore[:,:,t] + \n", " 1/M * Ht_ccc[:,:,t] + 1/M * Ht_comp_dcc[:,:,t] + 1/M* Ht_comp_ewma[:,:,t])\n", "\n", "min_final_w = minvar(H_mean_comb)\n", "min_final_eqw = minvar(H_min_eqw_comb)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "cut_ind = valid_wind\n", "consid_models=8\n", "m_weights= (min_final_w[cut_ind:],min_w_svech[cut_ind:], min_w_ore[cut_ind:],\n", " min_w_ccc[cut_ind:],min_w_dcc[cut_ind:], min_w_ewma[cut_ind:],equal_weights[cut_ind:],\n", " min_final_eqw[cut_ind:])\n", "m_names = ('COMBINATION','SVECH', 'ORE','CCC','DCC','EWMA','EQUAL WEIGHTS STOCKS', 'EQUAL WEIGHTS MODELS')\n", "\n", "see_perf(m_weights[:consid_models],m_names[:consid_models],ret[cut_ind:], False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## performance of both aproches combined" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cut_ind = valid_wind\n", "consid_models=20\n", "\n", "cut_start = valid_wind\n", "\n", "cut_end = T\n", "M_consid = 8\n", "m_weights= (equal_weights[cut_ind:cut_end],mean_w_svech[cut_start:cut_end], mean_w_ore[cut_start:cut_end],\n", " mean_w_ccc[cut_start:cut_end],mean_w_dcc[cut_start:cut_end], mean_w_ewma[cut_start:cut_end], \n", " min_w_svech[cut_start:cut_end], min_w_ore[cut_start:cut_end],\n", " min_w_ccc[cut_start:cut_end],min_w_dcc[cut_start:cut_end], min_w_ewma[cut_start:cut_end],\n", " )\n", "m_names = ('EQUAL WEIGHTS STOCKS','SVECH mean-var', 'ORE mean-var','CCC mean-var','DCC mean-var',\n", " 'EWMA mean-var','SVECH minvar', 'ORE minvar','CCC minvar','DCC minvar',\n", " 'EWMA minvar')\n", "\n", "candidate_perf = see_perf(m_weights[:consid_models],m_names[:consid_models],ret[cut_ind:cut_end], False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#candidate_perf.savefig('candgain.pdf')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Grid estimation\n", "used in numerical study" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Candidate models " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model_names = ('svech','ccc','dcc','ewma')\n", "M = len(model_names)\n", "\n", "Ht_svech= Ht_comp_svech\n", "Ht_ore = Ht_comp_ore\n", "Ht_dcc = Ht_comp_dcc\n", "Ht_ewma = Ht_comp_ewma\n", "\n", "\n", "H_eqw_comb = np.zeros((N,N,T)) #equal weights combination\n", "for t in range(T):\n", " H_eqw_comb[:,:,t] = (1/M* Ht_comp_svech[:,:,t] + \\\n", " 1/M * Ht_ccc[:,:,t] +1/M * Ht_comp_dcc[:,:,t] + 1/M* Ht_comp_ewma[:,:,t])\n", " \n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "perf_cand=[] \n", "validation_win=valid_wind\n", "deltas = np.ones(T)\n", "equal_weights = np.ones((T,N))*1/N\n", "delta_undisct =np.ones(T)\n", "delta = 1 \n", "\n", "eq_dperf, eq_ret, eq_ret_disc = port_return(equal_weights,ret,deltas)\n", "\n", "avg_turn = calc_turnover(equal_weights).mean()\n", "avg_dret = calc_cum_perf(eq_dperf[validation_win:])[-1]**(1/(T-validation_win))-1 \n", " \n", "undisc_var = np.var(eq_dperf[validation_win:])\n", "undisc_downvar = port_downside_variance(eq_dperf[validation_win:],eq_ret[validation_win:],\n", " delta_undisct,T-validation_win-1)\n", "undisc_beta = port_beta(eq_dperf[validation_win:], eq_ret[validation_win:],\n", " SP_ret[validation_win:],delta_undisct[validation_win:],T-validation_win-1)\n", " \n", "avg_dvol = np.sqrt(undisc_var)\n", "avg_beta =undisc_beta[-1]\n", " \n", "undisc_SR = sharpe_ratio(eq_ret[-1], undisc_var)\n", "undisc_SoR = sortino_ratio(eq_ret[-1], undisc_downvar[-1])\n", "undisc_TI = trainors_index(eq_ret[-1], undisc_beta[-1])\n", " \n", "avg_SR =undisc_SR\n", "avg_SoR = undisc_SoR\n", "avg_TI = undisc_TI \n", " \n", "\n", "\n", "\n", "perf_cand.append([\"-----\",\"Equal weights stocks\" ,avg_dret, avg_dvol, avg_turn,avg_beta, avg_SR,avg_SoR,avg_TI])\n", "\n", "for i in range(2):\n", " if i ==0:\n", " daily_exp_ret = 0.0008\n", " mu = rolling_mean(ret, 500)\n", " equal_weights_models = meanvar(H_eqw_comb, mu, daily_exp_ret)\n", " approach='meanvar'\n", " if i ==1:\n", " equal_weights_models = minvar(H_eqw_comb)\n", " approach='minvar'\n", " \n", " eq_dperf, eq_ret, eq_ret_disc = port_return(equal_weights_models,ret,deltas)\n", " #Save the metric we track\n", " eq_turn= calc_turnover(equal_weights_models)\n", " avg_turn = eq_turn[validation_win:].mean() \n", " avg_dret = calc_cum_perf(eq_dperf[validation_win:])[-1]**(1/(T-validation_win))-1 \n", " \n", " undisc_var = np.var(eq_dperf[validation_win:])\n", " undisc_downvar = port_downside_variance(eq_dperf[validation_win:],eq_ret[validation_win:],\n", " delta_undisct,T-validation_win-1)\n", " undisc_beta = port_beta(eq_dperf[validation_win:], eq_ret[validation_win:],\n", " SP_ret[validation_win:],delta_undisct[validation_win:],T-validation_win-1)\n", " \n", " avg_dvol = np.sqrt(undisc_var)\n", " avg_beta =undisc_beta[-1]\n", " \n", " undisc_SR = sharpe_ratio(eq_ret[-1], undisc_var)\n", " undisc_SoR = sortino_ratio(eq_ret[-1], undisc_downvar[-1])\n", " undisc_TI = trainors_index(eq_ret[-1], undisc_beta[-1])\n", " \n", " avg_SR =undisc_SR\n", " avg_SoR = undisc_SoR\n", " avg_TI = undisc_TI \n", " \n", " \n", " \n", " \n", " perf_cand.append([approach,\"Equal weights models\" ,avg_dret, avg_dvol, avg_turn,avg_beta, avg_SR,avg_SoR,avg_TI])\n", " \n", " \n", "perf_equal_m = pd.DataFrame(perf_cand,columns=['Approach',\"Model\" ,'avg_dret', 'avg_dvol', 'avg_turn','avg_beta', 'avg_SR','avg_SoR','avg_TI'])\n", "perf_equal_m[[\"avg_dret\",'avg_dvol','avg_TI']] = perf_equal_m[[\"avg_dret\",'avg_dvol','avg_TI']]*100\n", "perf_equal_m" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "# not dependent on delta or ni \n", "equal_weights = np.ones((T,N))*1/N\n", "validation_win = 300 # first observation to count in estimation\n", "deltas = np.ones(T)\n", "delta_undisct = deltas\n", "perf_cand=[]\n", "\n", "\n", "\n", "approaches = ['mean-var','minvar']\n", "model_names = ('svech','ore','ccc','dcc','ewma')\n", "cand_model_matrices = (Ht_comp_svech,Ht_comp_ore,Ht_ccc,Ht_comp_dcc,Ht_comp_ewma)\n", "#cand_model_matrices = (Ht_svech_svd,Ht_ore_svd,Ht_ccc,Ht_dcc_svd,Ht_ewma_svd)\n", "mu = rolling_mean(ret, 500) #column vector of expected returns\n", "\n", "all_w=[]\n", "\n", "\n", "for approach in approaches:\n", " for i,Ht in tqdm.tqdm(enumerate(cand_model_matrices)): \n", " if approach == 'mean-var':\n", " daily_exp_ret = 0.0008\n", " w = meanvar(Ht, mu, daily_exp_ret)\n", " else:\n", " w = minvar(Ht)\n", " \n", " all_w.append(w)\n", " turn = calc_turnover(w)\n", " dperf, cand_ret, ret_disc = port_return(w,ret,deltas) \n", " \n", " avg_turn = turn[validation_win:].mean() \n", " avg_dret = calc_cum_perf(dperf[validation_win:])[-1]**(1/(T-validation_win))-1 \n", "\n", " #Save the metric we track \n", " undisc_var = np.var(dperf[validation_win:])\n", " undisc_downvar = port_downside_variance(dperf[validation_win:],cand_ret[validation_win:],\n", " delta_undisct,T-validation_win-1)\n", " undisc_beta = port_beta(dperf[validation_win:], cand_ret[validation_win:],\n", " SP_ret[validation_win:],delta_undisct[validation_win:],T-validation_win-1)\n", " \n", " avg_dvol = np.sqrt(undisc_var)\n", " avg_beta =undisc_beta[-1]\n", " \n", " undisc_SR = sharpe_ratio(cand_ret[-1], undisc_var)\n", " undisc_SoR = sortino_ratio(cand_ret[-1], undisc_downvar[-1])\n", " undisc_TI = trainors_index(cand_ret[-1], undisc_beta[-1])\n", "\n", " \n", " avg_SR =undisc_SR\n", " avg_SoR = undisc_SoR\n", " avg_TI = undisc_TI \n", "\n", "\n", " perf_cand.append([approach,model_names[i], avg_dret, avg_dvol, avg_turn,avg_beta, avg_SR,avg_SoR,avg_TI])\n", "\n", "perf_cand_m = pd.DataFrame(perf_cand,columns=['Approach',\"Model\" ,'avg_dret', 'avg_dvol', 'avg_turn','avg_beta', 'avg_SR','avg_SoR','avg_TI'])\n", "perf_cand_m" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "perf_cand_m" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "perf_cand_m = pd.DataFrame(perf_cand,columns=['Approach',\"Model\" ,'avg_dret', 'avg_dvol', 'avg_turn','avg_beta', 'avg_SR','avg_SoR','avg_TI'])\n", "perf_cand_m[[\"avg_dret\",'avg_dvol','avg_TI']] = perf_cand_m[[\"avg_dret\",'avg_dvol','avg_TI']]*100\n", "perf_cand_m.to_csv('canidModelsf3mu8.csv',sep='&', float_format='%.3f',index=False)\n", "perf_cand_m" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "perf_cand_m" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Min-var " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "perf=[]\n", "validation_win = valid_wind # first observation to count in estimation\n", "mu = rolling_mean(ret, 500) #column vector of expected returns\n", "#approach = 'mean-var'\n", "approach = 'min-var'\n", "candm_weights_minvar=[]\n", "minvar_costadj_gain=[]\n", "\n", "delta_undisct =np.ones(T)\n", "\n", "metrics =['SR','SoR', 'TI']\n", "#metrics = ['SoR']\n", "delta_list =[1,0.99,0.95,0.9,0.8,0.7] \n", "ni_list=[1,2,5,10]\n", "\n", "if approach == 'mean-var':\n", " w_svech = meanvar(Ht_comp_svech, mu, daily_exp_ret) \n", " w_ore = meanvar(Ht_comp_ore, mu, daily_exp_ret)\n", " w_ccc= meanvar(Ht_ccc, mu, daily_exp_ret)\n", " w_dcc= meanvar(Ht_comp_dcc, mu, daily_exp_ret)\n", " w_ewma= meanvar(Ht_comp_ewma, mu, daily_exp_ret)\n", "else:\n", " w_svech = minvar(Ht_comp_svech)\n", " w_ore = minvar(Ht_comp_ore)\n", " w_ccc= minvar(Ht_ccc)\n", " w_dcc= minvar(Ht_comp_dcc)\n", " w_ewma= minvar(Ht_comp_ewma)\n", " \n", "\n", "for delta in tqdm.tqdm(delta_list):\n", " deltas = np.zeros(T)\n", " deltas[0] = 1\n", " deltas[1] = delta\n", " for t in range(2,T):\n", " deltas[t] = deltas[t-1]*delta\n", " \n", " for ni in ni_list:\n", " for metric in metrics: \n", "\n", " svech_dperf, svech_ret, svech_ret_disc = port_return(w_svech,ret,deltas) \n", " # ore_dperf, ore_ret, ore_ret_disc = port_return(w_ore,ret,deltas) #crossed out because of simillarity to ewma\n", " ccc_dperf, ccc_ret, ccc_ret_disc = port_return(w_ccc,ret,deltas)\n", " dcc_dperf, dcc_ret, dcc_ret_disc = port_return(w_dcc,ret,deltas)\n", " ewma_dperf, ewma_ret, ewma_ret_disc = port_return(w_ewma,ret,deltas)\n", " \n", "\n", " svech_var = port_variance(svech_dperf,svech_ret,deltas)\n", " # ore_var = port_variance(ore_dperf,ore_ret,deltas)\n", " dcc_var = port_variance(dcc_dperf,dcc_ret,deltas)\n", " ccc_var = port_variance(ccc_dperf,ccc_ret,deltas)\n", " ewma_var = port_variance(ewma_dperf,ewma_ret,deltas)\n", " \n", "\n", " svech_downvar = port_downside_variance(svech_dperf,svech_ret,deltas)\n", " # ore_downvar = port_downside_variance(ore_dperf,ore_ret,deltas)\n", " dcc_downvar = port_downside_variance(dcc_dperf,dcc_ret,deltas)\n", " ccc_downvar = port_downside_variance(ccc_dperf,ccc_ret,deltas)\n", " ewma_downvar = port_downside_variance(ewma_dperf,ewma_ret,deltas)\n", " \n", " \n", " svech_beta = port_beta(svech_dperf, svech_ret,SP_ret,deltas)\n", " # ore_beta = port_beta(ore_dperf, ore_ret,SP_ret,deltas)\n", " ccc_beta = port_beta(ccc_dperf, ccc_ret,SP_ret,deltas)\n", " dcc_beta = port_beta(dcc_dperf, dcc_ret,SP_ret,deltas)\n", " ewma_beta = port_beta(ewma_dperf, ewma_ret,SP_ret,deltas)\n", " \n", "\n", " svech_SR = sharpe_ratio(svech_ret_disc, svech_var)\n", " # ore_SR = sharpe_ratio(ore_ret_disc,ore_var)\n", " ccc_SR = sharpe_ratio(ccc_ret_disc,ccc_var)\n", " dcc_SR = sharpe_ratio(dcc_ret_disc,dcc_var)\n", " ewma_SR = sharpe_ratio(ewma_ret_disc, ewma_var)\n", " \n", "\n", " svech_SoR = sortino_ratio(svech_ret_disc, svech_downvar)\n", " # ore_SoR = sortino_ratio(ore_ret_disc,ore_downvar)\n", " ccc_SoR = sortino_ratio(ccc_ret_disc,ccc_downvar)\n", " dcc_SoR = sortino_ratio(dcc_ret_disc,dcc_downvar)\n", " ewma_SoR = sortino_ratio(ewma_ret_disc, ewma_downvar)\n", " \n", "\n", " \n", " svech_TI = trainors_index(svech_ret_disc, svech_beta)\n", " # ore_TI = trainors_index(ore_ret_disc, ore_beta)\n", " ccc_TI = trainors_index(ccc_ret_disc, ccc_beta)\n", " dcc_TI = trainors_index(dcc_ret_disc, dcc_beta)\n", " ewma_TI = trainors_index(ewma_ret_disc, ewma_beta)\n", " \n", "\n", " if metric =='SR':\n", " models = (svech_SR,ccc_SR,dcc_SR,ewma_SR)\n", " elif metric =='SoR':\n", " models = (svech_SoR,ccc_SoR,dcc_SoR,ewma_SoR)\n", " else: #TI chosen\n", " models = (svech_TI,ccc_TI,dcc_TI,ewma_TI)\n", "\n", " model_names =('svech','ccc','dcc','ewma')\n", " M = len(models)\n", " init_positions = np.ones(M) * 1/M\n", "\n", " stack = np.stack(arrays = models)\n", " stack= np.where(stack<0, 0, stack)\n", " stack = stack **ni\n", "\n", " cum_weights = stack.sum(axis = 0 )\n", " model_weights = stack / cum_weights\n", " model_weights[:,0] = init_positions\n", "\n", " indx= np.where(cum_weights<=0) #cases when all SRs are below zero - choose equal weights\n", " model_weights[:,indx] = 1/M\n", " model_weights[:,[0,1]] = 1/m # initial setup\n", "\n", " shifted_weights = np.roll(model_weights, -1, axis=1)\n", "\n", " shifted_weights[:,0] = init_positions\n", "\n", " shifted_weights = np.insert(model_weights, 0, init_positions, axis=1) # initial setup\n", " shifted_weights = np.delete(shifted_weights,obj = -1, axis =1)\n", "\n", " H_comb = np.zeros((N,N,T))\n", " \n", " M=shifted_weights.shape[0]\n", "\n", " for t in range(T):\n", " H_comb[:,:,t] = (shifted_weights[0,t]* Ht_comp_svech[:,:,t] + \\\n", " shifted_weights[1,t] * Ht_ccc[:,:,t]+shifted_weights[2,t] * Ht_comp_dcc[:,:,t] + \\\n", " shifted_weights[3,t]* Ht_comp_ewma[:,:,t])\n", " if approach == 'mean-var':\n", " daily_exp_ret = 0.0005\n", " final_w = meanvar(H_comb, mu, daily_exp_ret)\n", " else:\n", " final_w = minvar(H_comb)\n", " \n", " \n", " final_turn = calc_turnover(final_w)\n", " final_dperf, final_ret, final_ret_disc = port_return(final_w,ret,deltas) \n", " final_gain_costadj = port_return_costadj(final_w[validation_win:],ret[validation_win:],turn[validation_win:])\n", " minvar_costadj_gain.append(final_gain_costadj)\n", " avg_turn = final_turn[validation_win:].mean() \n", " avg_dret = calc_cum_perf(final_dperf[validation_win:])[-1]**(1/(T-validation_win))-1 \n", "\n", " #Save the metric we track \n", " undisc_var = np.var(final_dperf[validation_win:])\n", " undisc_downvar = port_downside_variance(final_dperf[validation_win:],final_ret[validation_win:],\n", " delta_undisct,T-validation_win-1)\n", " undisc_beta = port_beta(final_dperf[validation_win:], final_ret[validation_win:],\n", " SP_ret[validation_win:],delta_undisct[validation_win:],T-validation_win-1)\n", " \n", " avg_dvol = np.sqrt(undisc_var)\n", " avg_beta =undisc_beta[-1]\n", " \n", " undisc_SR = sharpe_ratio(final_ret[-1], undisc_var)\n", " undisc_SoR = sortino_ratio(final_ret[-1], undisc_downvar[-1])\n", " undisc_TI = trainors_index(final_ret[-1], undisc_beta[-1])\n", " \n", " avg_SR =undisc_SR\n", " avg_SoR = undisc_SoR\n", " avg_TI = undisc_TI \n", " \n", " perf.append([approach,metric,delta, ni, avg_dret, avg_dvol, avg_turn,avg_beta, avg_SR,avg_SoR,avg_TI])\n", " candm_weights_minvar.append(shifted_weights)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "minvar_df = pd.DataFrame(perf,columns=['approach','metric','delta', 'ni', 'avg_dret', 'avg_dvol', 'avg_turn','avg_beta', 'avg_SR','avg_SoR','avg_TI'])\n", "minvar_df" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "minvar_df = pd.DataFrame(perf,columns=['approach','metric','delta', 'ni', 'avg_dret', 'avg_dvol', 'avg_turn','avg_beta', 'avg_SR','avg_SoR','avg_TI'])\n", "minvar_df['new']= minvar_df['approach'] + minvar_df['delta'].apply(lambda x: \"(\"+str(x))+minvar_df['ni'].apply(lambda x: \" ,\"+str(x))+ minvar_df['metric'].apply(lambda x: \" ,\"+str(x)+\")\")\n", "minvar_output = minvar_df.drop(labels = ['approach','metric','delta','ni'],axis= 1 )\n", "cols =minvar_output.columns.tolist()\n", "cols = cols[-1:] + cols[:-1]\n", "minvar_output[[\"avg_dret\",'avg_dvol','avg_TI']] = minvar_output[[\"avg_dret\",'avg_dvol','avg_TI']]*100\n", "minvar_output = minvar_output[cols]\n", "minvar_output" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "minvar_output.sort_values(by='avg_dret', ascending= False).iloc[:10].to_csv('TOP10minvarModelsf3.csv',sep='&', float_format='%.3f',index=False)\n", "minvar_output.sort_values(by='avg_dret', ascending= False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#minvar_output.to_csv('minvarModelsf3.csv',sep='&', float_format='%.3f',index=False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "a=[]\n", "for elem in minvar_costadj_gain:\n", " a.append(elem[4])\n", "plt.plot(a)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mean-var" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "perf=[]\n", "candm_weights_meanvar=[]\n", "meanvar_costadj_gain=[]\n", "\n", "delta_undisct =np.ones(T)\n", "\n", "validation_win = valid_wind # first observation to count in estimation\n", "mu = rolling_mean(ret, 500) #column vector of expected returns\n", "approach = 'mean-var'\n", "\n", "metrics =['SR','SoR', 'TI']\n", "#metrics = ['SoR']\n", "delta_list =[1,0.99,0.95,0.9,0.8,0.7] \n", "ni_list=[1,2,5,10]\n", "#delta_list =[0.7,0.5,0.2]\n", "if approach == 'mean-var':\n", " w_svech = meanvar(Ht_comp_svech, mu, daily_exp_ret) \n", " #w_ore = meanvar(Ht_comp_ore, mu, daily_exp_ret)\n", " w_ccc= meanvar(Ht_ccc, mu, daily_exp_ret)\n", " w_dcc= meanvar(Ht_comp_dcc, mu, daily_exp_ret)\n", " w_ewma= meanvar(Ht_comp_ewma, mu, daily_exp_ret)\n", "else:\n", " w_svech = minvar(Ht_comp_svech)\n", " #w_ore = minvar(Ht_comp_ore)\n", " w_ccc= minvar(Ht_ccc)\n", " w_dcc= minvar(Ht_comp_dcc)\n", " w_ewma= minvar(Ht_comp_ewma)\n", " \n", "\n", "for delta in tqdm.tqdm(delta_list):\n", " deltas = np.zeros(T)\n", " deltas[0] = 1\n", " deltas[1] = delta\n", " for t in range(2,T):\n", " deltas[t] = deltas[t-1]*delta\n", " \n", " for ni in ni_list:\n", " for metric in metrics: \n", "\n", "\n", " #calc model metrics \n", " svech_dperf, svech_ret, svech_ret_disc = port_return(w_svech,ret,deltas) \n", " # ore_dperf, ore_ret, ore_ret_disc = port_return(w_ore,ret,deltas)\n", " ccc_dperf, ccc_ret, ccc_ret_disc = port_return(w_ccc,ret,deltas)\n", " dcc_dperf, dcc_ret, dcc_ret_disc = port_return(w_dcc,ret,deltas)\n", " ewma_dperf, ewma_ret, ewma_ret_disc = port_return(w_ewma,ret,deltas)\n", " \n", "\n", " svech_var = port_variance(svech_dperf,svech_ret,deltas)\n", " # ore_var = port_variance(ore_dperf,ore_ret,deltas)\n", " dcc_var = port_variance(dcc_dperf,dcc_ret,deltas)\n", " ccc_var = port_variance(ccc_dperf,ccc_ret,deltas)\n", " ewma_var = port_variance(ewma_dperf,ewma_ret,deltas)\n", " \n", "\n", " svech_downvar = port_downside_variance(svech_dperf,svech_ret,deltas)\n", " # ore_downvar = port_downside_variance(ore_dperf,ore_ret,deltas)\n", " dcc_downvar = port_downside_variance(dcc_dperf,dcc_ret,deltas)\n", " ccc_downvar = port_downside_variance(ccc_dperf,ccc_ret,deltas)\n", " ewma_downvar = port_downside_variance(ewma_dperf,ewma_ret,deltas)\n", " \n", " \n", " svech_beta = port_beta(svech_dperf, svech_ret,SP_ret,deltas)\n", " # ore_beta = port_beta(ore_dperf, ore_ret,SP_ret,deltas)\n", " ccc_beta = port_beta(ccc_dperf, ccc_ret,SP_ret,deltas)\n", " dcc_beta = port_beta(dcc_dperf, dcc_ret,SP_ret,deltas)\n", " ewma_beta = port_beta(ewma_dperf, ewma_ret,SP_ret,deltas)\n", " \n", "\n", " svech_SR = sharpe_ratio(svech_ret_disc, svech_var)\n", " # ore_SR = sharpe_ratio(ore_ret_disc,ore_var)\n", " ccc_SR = sharpe_ratio(ccc_ret_disc,ccc_var)\n", " dcc_SR = sharpe_ratio(dcc_ret_disc,dcc_var)\n", " ewma_SR = sharpe_ratio(ewma_ret_disc, ewma_var)\n", " \n", "\n", " svech_SoR = sortino_ratio(svech_ret_disc, svech_downvar)\n", " # ore_SoR = sortino_ratio(ore_ret_disc,ore_downvar)\n", " ccc_SoR = sortino_ratio(ccc_ret_disc,ccc_downvar)\n", " dcc_SoR = sortino_ratio(dcc_ret_disc,dcc_downvar)\n", " ewma_SoR = sortino_ratio(ewma_ret_disc, ewma_downvar)\n", " \n", "\n", " \n", " svech_TI = trainors_index(svech_ret_disc, svech_beta)\n", " #ore_TI = trainors_index(ore_ret_disc, ore_beta)\n", " ccc_TI = trainors_index(ccc_ret_disc, ccc_beta)\n", " dcc_TI = trainors_index(dcc_ret_disc, dcc_beta)\n", " ewma_TI = trainors_index(ewma_ret_disc, ewma_beta)\n", " \n", "\n", " if metric =='SR':\n", " models = (svech_SR,ccc_SR,dcc_SR,ewma_SR)\n", " elif metric =='SoR':\n", " models = (svech_SoR,ccc_SoR,dcc_SoR,ewma_SoR)\n", " else: #TI chosen\n", " models = (svech_TI,ccc_TI,dcc_TI,ewma_TI)\n", "\n", " model_names =('svech','ccc','dcc','ewma')\n", " M = len(models)\n", " init_positions = np.ones(M) * 1/M\n", "\n", " stack = np.stack(arrays = models)\n", " stack= np.where(stack<0, 0, stack)\n", " stack = stack **ni\n", "\n", " cum_weights = stack.sum(axis = 0 )\n", " model_weights = stack / cum_weights\n", " model_weights[:,0] = init_positions\n", "\n", " indx= np.where(cum_weights<=0) #cases when all SRs are below zero - choose equal weights\n", " model_weights[:,indx] = 1/M\n", " model_weights[:,[0,1]] = 1/len(models) # initial setup\n", "\n", " shifted_weights = np.roll(model_weights, -1, axis=1)\n", "\n", " shifted_weights[:,0] = init_positions\n", "\n", " shifted_weights = np.insert(model_weights, 0, init_positions, axis=1) # initial setup\n", " shifted_weights = np.delete(shifted_weights,obj = -1, axis =1)\n", "\n", " H_comb = np.zeros((N,N,T))\n", " \n", " M=shifted_weights.shape[0]\n", "\n", " for t in range(T):\n", " H_comb[:,:,t] = (shifted_weights[0,t]* Ht_comp_svech[:,:,t] + \\\n", " shifted_weights[1,t] * Ht_ccc[:,:,t]+shifted_weights[2,t] * Ht_comp_dcc[:,:,t] + \\\n", " shifted_weights[3,t]* Ht_comp_ewma[:,:,t])\n", " if approach == 'mean-var':\n", " daily_exp_ret = 0.0008\n", " final_w = meanvar(H_comb, mu, daily_exp_ret)\n", " else:\n", " final_w = minvar(H_comb)\n", " \n", " \n", " final_turn = calc_turnover(final_w)\n", " final_dperf, final_ret, final_ret_disc = port_return(final_w,ret,deltas) \n", " \n", " avg_turn = final_turn[validation_win:].mean() \n", " avg_dret = calc_cum_perf(final_dperf[validation_win:])[-1]**(1/(T-validation_win))-1 \n", "\n", " #Save the metric we track\n", " \n", " \n", " undisc_var = np.var(final_dperf[validation_win:])\n", " undisc_downvar = port_downside_variance(final_dperf[validation_win:],final_ret[validation_win:],\n", " delta_undisct,T-validation_win-1)\n", " undisc_beta = port_beta(final_dperf[validation_win:], final_ret[validation_win:],\n", " SP_ret[validation_win:],delta_undisct[validation_win:],T-validation_win-1)\n", " \n", " final_gain_costadj = port_return_costadj(final_w[validation_win:],ret[validation_win:],turn[validation_win:])\n", " meanvar_costadj_gain.append(final_gain_costadj)\n", " \n", " avg_dvol = np.sqrt(undisc_var)\n", " avg_beta =undisc_beta[-1]\n", " \n", " undisc_SR = sharpe_ratio(final_ret[-1], undisc_var)\n", " undisc_SoR = sortino_ratio(final_ret[-1], undisc_downvar[-1])\n", " undisc_TI = trainors_index(final_ret[-1], undisc_beta[-1])\n", " \n", " avg_SR =undisc_SR\n", " avg_SoR = undisc_SoR\n", " avg_TI = undisc_TI\n", " perf.append([approach,metric,delta, ni, avg_dret, avg_dvol, avg_turn,avg_beta, avg_SR,avg_SoR,avg_TI])\n", " candm_weights_meanvar.append(shifted_weights)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "meanvar_df = pd.DataFrame(perf,columns=['approach','metric','delta', 'ni', 'avg_dret', 'avg_dvol', 'avg_turn','avg_beta', 'avg_SR','avg_SoR','avg_TI'])\n", "meanvar_df['new']= meanvar_df['approach'] + meanvar_df['delta'].apply(lambda x: \"(\"+str(x))+meanvar_df['ni'].apply(lambda x: \" ,\"+str(x))+ meanvar_df['metric'].apply(lambda x: \" ,\"+str(x)+\")\")\n", "meanvar_output = meanvar_df.drop(labels = ['approach','metric','delta','ni'],axis= 1 )\n", "cols = meanvar_output.columns.tolist()\n", "cols = cols[-1:] + cols[:-1]\n", "meanvar_output[[\"avg_dret\",'avg_dvol','avg_TI']] = meanvar_output[[\"avg_dret\",'avg_dvol','avg_TI']]*100\n", "meanvar_output = meanvar_output[cols]\n", "meanvar_output" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "meanvar_output.sort_values(by='avg_dret',ascending=False).iloc[:10].to_csv('meanvarModelsf3Top10.csv',sep='&', float_format='%.3f',index=False)\n", "meanvar_output.sort_values(by='avg_dret',ascending=False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#meanvar_output.to_csv('meanvarModelsf3.csv',sep='&', float_format='%.3f',index=False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "a= pd.DataFrame(meanvar_costadj_gain, columns=[0,5,10,15,20], index = meanvar_output['new'])\n", "a= a.sort_values(by=0, ascending= False)\n", "a.to_csv('meancumul.csv',sep='&', float_format='%.3f',index=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "a= pd.DataFrame(minvar_costadj_gain, columns=[0,5,10,15,20], index = minvar_output['new'])\n", "a= a.sort_values(by=0, ascending= False)\n", "a.to_csv('mincumul.csv',sep='&', float_format='%.3f',index=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pd.options.display.float_format = '{:,.3f}'.format\n", "a= pd.DataFrame(minvar_costadj_gain, columns=[0,5,10,15,20], index = minvar_output['new'])\n", "\n", "a" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "meanvar_output" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [], "scrolled": false }, "outputs": [], "source": [ "# figure showed in thesis\n", "model_names=('BEKK', 'CCC', 'DCC', 'EWMA')\n", "fig,axe = plt.subplots(2,3,figsize = (22,10))\n", "x_pos = np.arange(len(model_names))\n", "lst =[1,2,3,21,22,23]\n", "for i,g in enumerate(lst): # range(len(candm_weights_meanvar)):\n", " wei = pd.DataFrame(candm_weights_minvar[g].T)\n", " end = T\n", " s=300\n", " x=range(0,5) \n", " data= [wei.iloc[s:end,0],wei.iloc[s:end,1], wei.iloc[s:end,2], wei.iloc[s:end,3]]\n", " #plt.boxplot(wei.iloc[:end,0],wei.iloc[:end,1], wei.iloc[:end,2], wei.iloc[:end,3],wei.iloc[:end,4])\n", " \n", " axe[i //3,i %3].boxplot(data , vert = False, showfliers=False) \n", " axe[i //3,i %3].set_yticklabels( model_names[:4])\n", " axe[i //3,i %3].set_xlim(-0.05,1.05) \n", " if i%3 == 0:\n", " riskmetric ='SR'\n", " elif i%3 == 1:\n", " riskmetric ='SoR'\n", " else: riskmetric ='TI'\n", " title = riskmetric+ str(r\" $\\delta=$\" + str(delta_list[g//12])) + r\" $\\eta=$\"+ str(ni_list[((g//6)%4)])\n", " axe[i //3,i %3].set_title(title)\n", " # axe[g //3,g %3].boxplot(data , vert = False, showfliers=False) \n", " #axe[g //3,g %3].set_yticklabels( model_names[:4])\n", " #axe[g //3,g %3].set_xlim(-0.05,1.05) \n", " #title = str(r\" $\\delta=$\" + str(delta_list[g//9])) + r\" $\\eta=$\"+ str(ni_list[((g//3)%3)])\n", " #axe[g //3,g %3].set_title(title)\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#fig.savefig('minvarweights.pdf')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "# figure showed in thesis\n", "model_names=('BEKK', 'CCC', 'DCC', 'EWMA')\n", "fig,axe = plt.subplots(4,3,figsize = (30,35))\n", "x_pos = np.arange(len(model_names))\n", "lst =range(12,24)\n", "for i,g in enumerate(lst): # range(len(candm_weights_meanvar)):\n", " wei = pd.DataFrame(candm_weights_meanvar[g].T)\n", " end = T\n", " s=300\n", " x=range(0,5) \n", " data= [wei.iloc[s:end,0],wei.iloc[s:end,1], wei.iloc[s:end,2], wei.iloc[s:end,3]]\n", " #plt.boxplot(wei.iloc[:end,0],wei.iloc[:end,1], wei.iloc[:end,2], wei.iloc[:end,3],wei.iloc[:end,4])\n", " \n", " axe[i //3,i %3].boxplot(data , vert = False, showfliers=False) \n", " axe[i //3,i %3].set_yticklabels( model_names[:4])\n", " axe[i //3,i %3].set_xlim(-0.05,1.05) \n", " if i%3 == 0:\n", " riskmetric ='SR'\n", " elif i%3 == 1:\n", " riskmetric ='SoR'\n", " else: riskmetric ='TI'\n", " title = riskmetric+ str(r\" $\\delta=$\" + str(delta_list[g//12])) + r\" $\\eta=$\"+ str(ni_list[((g//3)%4)])\n", " axe[i //3,i %3].set_title(title, fontsize = 40)\n", " # axe[g //3,g %3].boxplot(data , vert = False, showfliers=False) \n", " #axe[g //3,g %3].set_yticklabels( model_names[:4])\n", " #axe[g //3,g %3].set_xlim(-0.05,1.05) \n", " #title = str(r\" $\\delta=$\" + str(delta_list[g//9])) + r\" $\\eta=$\"+ str(ni_list[((g//3)%3)])\n", " #axe[g //3,g %3].set_title(title)\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ " #fig.savefig('meanvarweightsfull.pdf')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "# weight of all portfolios - not showed in thesis because of the size\n", "model_names=('BEKK', 'CCC', 'DCC', 'EWMA')\n", "fig,axe = plt.subplots(int(72/3),3,figsize = (20,90))\n", "x_pos =[ np.arange(len(model_names))]\n", "lst =range(0,72)\n", "for i,g in enumerate(lst): # range(len(candm_weights_meanvar)):\n", " wei = pd.DataFrame(candm_weights_meanvar[g].T)\n", " end = T\n", " s=300\n", " x=range(0,5) \n", " data= [wei.iloc[s:end,0],wei.iloc[s:end,1], wei.iloc[s:end,2], wei.iloc[s:end,3]]\n", " #plt.boxplot(wei.iloc[:end,0],wei.iloc[:end,1], wei.iloc[:end,2], wei.iloc[:end,3],wei.iloc[:end,4])\n", " \n", " axe[i //3,i %3].boxplot(data , vert = False, showfliers=False) \n", " axe[i //3,i %3].set_yticklabels( model_names[:4])\n", " axe[i //3,i %3].set_xlim(-0.05,1.05) \n", " if i%3 == 0:\n", " riskmetric ='SR'\n", " elif i%3 == 1:\n", " riskmetric ='SoR'\n", " else: riskmetric ='TI'\n", " title = riskmetric+ str(r\" $\\delta=$\" + str(delta_list[g//12])) + r\" $\\eta=$\"+ str(ni_list[((g//3)%4)])\n", " axe[i //3,i %3].set_title(title, fontsize = 10)\n", " # axe[g //3,g %3].boxplot(data , vert = False, showfliers=False) \n", " #axe[g //3,g %3].set_yticklabels( model_names[:4])\n", " #axe[g //3,g %3].set_xlim(-0.05,1.05) \n", " #title = str(r\" $\\delta=$\" + str(delta_list[g//9])) + r\" $\\eta=$\"+ str(ni_list[((g//3)%3)])\n", " #axe[g //3,g %3].set_title(title)\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# fig.savefig('minvarweightsfull.pdf')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Testing the model specification\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "estimates_names =[\n", "'EWMA_est',\n", "'EWMA_svd_est',\n", "'EWMA_comp_est',\n", "'EWMA_comp_est_full',\n", "'CCC_est',\n", "'DCC_est',\n", "'DCC_svd_est',\n", "'DCC_comp_est',\n", "'DCC_comp_est_full',\n", "'SVECH_est',\n", "'SVECH_svd_est',\n", "'SVECH_comp_est',\n", "'SVECH_comp_est_full',\n", "'ORE_est',\n", "'ORE_svd_est',\n", "'ORE_comp_est',\n", "'ORE_comp_est_full']\n", "\n", "\n", "matrices = [Ht_ewma, \n", "Ht_ewma_svd,\n", "Ht_comp_ewma,\n", "Ht_comp_ewma_full,\n", "Ht_ccc,\n", "Ht_dcc,\n", "Ht_dcc_svd,\n", "Ht_comp_dcc,\n", "Ht_comp_dcc_full,\n", "Ht_svech,\n", "Ht_svech_svd,\n", "Ht_comp_svech,\n", "Ht_comp_svech_full,\n", "Ht_ore,\n", "Ht_ore_svd,\n", "Ht_comp_ore,\n", "Ht_comp_ore_full]\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def get_trace(C0_inv,Ci):\n", " mat = np.matmul(Ci.T,C0_inv)\n", " mat = np.matmul(mat, Ci)\n", " mat=np.matmul(mat, C0_inv)\n", " \n", " return np.trace(mat)\n", " \n", "\n", "def Port_stat(H,ret,h):\n", " T= H.shape[2]\n", " N= H.shape[0]\n", " H_chol_inv = np.zeros((N,N,T))\n", " H_inv_ret = np.zeros((N,N,T))\n", " C=np.zeros((N,N,h+1))\n", " trace = np.zeros(h+1)\n", " \n", " test_stat =0\n", " \n", " for t in range(T):\n", " H_chol_inv[:,:,t]= np.linalg.inv(np.linalg.cholesky(H[:,:,t]))\n", " H_inv_ret[:,:,t] =np.array(np.matmul(H_chol_inv[:,:,t], ret[t,:]),ndmin=2).T\n", " for j in range(h+1):\n", " for t in range(j,T): \n", " C[:,:,j]=C[:,:,j] + np.matmul(H_inv_ret[:,:,t],H_inv_ret[:,:,t-j].T)\n", " \n", " C[:,:,j] = 1/T*C[:,:,j]\n", " \n", " C0_inv = np.linalg.inv(C[:,:,0])\n", " \n", " for j in range(1,h+1):\n", " trace[j] = get_trace(C0_inv, C[:,:,j])\n", " test_stat =test_stat+ 1/(T-j)*trace[j]\n", " \n", " return T**2*test_stat" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def Portmantau_stat2(H,ret,h):\n", " T= H.shape[2]\n", " N= H.shape[0]\n", " H_inv = np.zeros((N,N,T))\n", " H_inv_ret = np.zeros((N,N,T))\n", " rhr = np.zeros(T)\n", " test_stat= np.zeros(h+1)\n", " \n", " for t in range(T): \n", " H_inv[:,:,t] =np.linalg.inv(H[:,:,t])\n", " rhr[t] = np.matmul(np.matmul(np.array(ret[t,:],ndmin=2), H_inv[:,:,t]),np.array(ret[t,:],ndmin=2).T) \n", " test_statj=0\n", " denomj=0\n", " \n", " for j in range(1,h+1): \n", " for t in range(j+1,T+1):\n", " test_statj += ((rhr[t-1]-N)*(rhr[t-j-1])-N)\n", " denomj += (rhr[t-1]-N)**2 \n", " test_stat[j]= test_statj/denomj \n", " return (test_stat**2).sum()*T" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "port_tests1=[]\n", "\n", "for matrix in matrices:\n", " port_tests1.append(Port_stat(matrix[:,:,300:],ret[300:,:],10))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for i,elem in enumerate(port_tests1):\n", " print(elem,estimates_names[i])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "port_tests2=[]\n", "for matrix in matrices:\n", " port_tests2.append(Portmantau_stat2(matrix[:,:,300:],ret[300:,:],10))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for i,elem in enumerate(port_tests2):\n", " print(elem,estimates_names[i])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "port_tests=[]\n", "\n", "matrices = [Ht_ewma, \n", "Ht_ewma_svd,\n", "#Ht_comp_ewma,\n", "Ht_ccc,\n", "Ht_dcc,\n", "Ht_dcc_svd,\n", "#Ht_comp_dcc,\n", "\n", "Ht_svech,\n", "Ht_svech_svd,\n", "#Ht_comp_svech,\n", "\n", "Ht_ore,\n", "Ht_ore_svd]\n", "#Ht_comp_ore]\n", "\n", "\n", "for matrix in matrices:\n", " port_tests.append(Portmantau_stat2(matrix[:,:,300:],ret[300:,:],10))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for elem in port_tests:\n", " print(elem)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": { "height": "calc(100% - 180px)", "left": "10px", "top": "150px", "width": "299.2px" }, "toc_section_display": true, "toc_window_display": true }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }