[[
wikihub
]]
Search
⌘K
Explore
Activity
People
For Agents
Sign in
Explore
Activity
People
For Agents
Sign in
×
@jemoka / Jemoka Knowledge Base / raw/course/nus_econ320/kbhnus_econ320_linearity_tests.md
Suggest edit
Cancel
Submit suggestion
Title
Name
Note
--- title: "NUS-ECON320 Linearity Tests" source: https://www.jemoka.com/posts/kbhnus_econ320_linearity_tests/ --- Let’s begin. We want to create test for the linearity of a few assets, for whether or not they follow the CAPM. Note that we will be using the Sharpe-Linter version of CAPM: \begin{equation} E[R_{i}-R_{f}] = \beta_{im} E[(R_{m}-R_{f})] \end{equation} \begin{equation} \beta_{im} := \frac{Cov[(R_{i}-R_{f}),(R_{m}-R_{f})]}{Var[R_{m}-R_{f}]} \end{equation} Recall that we declare \(R_{f}\) (the risk-free rate) to be non-stochastic. Let us begin. We will create a generic function to analyze some given stock. We will first import our utilities import pandas as pd import numpy as np Let’s load the data from our market (NYSE) as well as our 10 year t-bill data. t_bill = pd.read_csv("./linearity_test_data/10yrTBill.csv") nyse = pd.read_csv("./linearity_test_data/NYSEComposite.csv") nyse.head() Date Close 0 11/7/2013 16:00:00 9924.37 1 11/8/2013 16:00:00 10032.14 2 11/11/2013 16:00:00 10042.95 3 11/12/2013 16:00:00 10009.84 4 11/13/2013 16:00:00 10079.89 Excellent. Let’s load in the data for that stock. def load_stock(stock): return pd.read_csv(f"./linearity_test_data/{stock}.csv") load_stock("LMT").head() Date Close 0 11/7/2013 16:00:00 136.20 1 11/8/2013 16:00:00 138.11 2 11/11/2013 16:00:00 137.15 3 11/12/2013 16:00:00 137.23 4 11/13/2013 16:00:00 137.26 And now, let’s load all three stocks, then concatenate them all into a big-ol DataFrame. # load data df = { "Date": nyse.Date, "NYSE": nyse.Close, "TBill": t_bill.Close, "LMT": load_stock("LMT").Close, "TWTR": load_stock("TWTR").Close, "MCD": load_stock("MCD").Close } # convert to dataframe df = pd.DataFrame(df) # drop empty df.dropna(inplace=True) df Date NYSE TBill LMT TWTR MCD 0 11/7/2013 16:00:00 9924.37 26.13 136.20 44.90 97.20 1 11/8/2013 16:00:00 10032.14 27.46 138.11 41.65 97.01 2 11/11/2013 16:00:00 10042.95 27.51 137.15 42.90 97.09 3 11/12/2013 16:00:00 10009.84 27.68 137.23 41.90 97.66 4 11/13/2013 16:00:00 10079.89 27.25 137.26 42.60 98.11 ... ... ... ... ... ... ... 2154 10/24/2022 16:00:00 14226.11 30.44 440.11 39.60 252.21 2155 10/25/2022 16:00:00 14440.69 31.56 439.30 39.30 249.28 2156 10/26/2022 16:00:00 14531.69 33.66 441.11 39.91 250.38 2157 10/27/2022 16:00:00 14569.90 34.83 442.69 40.16 248.36 2158 10/28/2022 16:00:00 14795.63 33.95 443.41 39.56 248.07 [2159 rows x 6 columns] Excellent. Now, let’s convert all of these values into daily log-returns (we don’t really care about the actual pricing.) log_returns = df[["NYSE", "TBill", "LMT", "TWTR", "MCD"]].apply(np.log, inplace=True) df.loc[:, ["NYSE", "TBill", "LMT", "TWTR", "MCD"]] = log_returns df Date NYSE TBill LMT TWTR MCD 0 11/7/2013 16:00:00 9.202749 3.263084 4.914124 3.804438 4.576771 1 11/8/2013 16:00:00 9.213549 3.312730 4.928050 3.729301 4.574814 2 11/11/2013 16:00:00 9.214626 3.314550 4.921075 3.758872 4.575638 3 11/12/2013 16:00:00 9.211324 3.320710 4.921658 3.735286 4.581492 4 11/13/2013 16:00:00 9.218298 3.305054 4.921877 3.751854 4.586089 ... ... ... ... ... ... ... 2154 10/24/2022 16:00:00 9.562834 3.415758 6.087025 3.678829 5.530262 2155 10/25/2022 16:00:00 9.577805 3.451890 6.085183 3.671225 5.518577 2156 10/26/2022 16:00:00 9.584087 3.516310 6.089294 3.686627 5.522980 2157 10/27/2022 16:00:00 9.586713 3.550479 6.092870 3.692871 5.514879 2158 10/28/2022 16:00:00 9.602087 3.524889 6.094495 3.677819 5.513711 [2159 rows x 6 columns] We will now calculate the log daily returns. But before—the dates are no longer relavent here so we drop them. df Date NYSE TBill LMT TWTR MCD 0 11/7/2013 16:00:00 9.202749 3.263084 4.914124 3.804438 4.576771 1 11/8/2013 16:00:00 9.213549 3.312730 4.928050 3.729301 4.574814 2 11/11/2013 16:00:00 9.214626 3.314550 4.921075 3.758872 4.575638 3 11/12/2013 16:00:00 9.211324 3.320710 4.921658 3.735286 4.581492 4 11/13/2013 16:00:00 9.218298 3.305054 4.921877 3.751854 4.586089 ... ... ... ... ... ... ... 2154 10/24/2022 16:00:00 9.562834 3.415758 6.087025 3.678829 5.530262 2155 10/25/2022 16:00:00 9.577805 3.451890 6.085183 3.671225 5.518577 2156 10/26/2022 16:00:00 9.584087 3.516310 6.089294 3.686627 5.522980 2157 10/27/2022 16:00:00 9.586713 3.550479 6.092870 3.692871 5.514879 2158 10/28/2022 16:00:00 9.602087 3.524889 6.094495 3.677819 5.513711 [2159 rows x 6 columns] And now, the log returns! We will shift this data by one column and subtract. returns = df.drop(columns=["Date"]) - df.drop(columns=["Date"]).shift(1) returns.dropna(inplace=True) returns NYSE TBill LMT TWTR MCD 1 0.010801 0.049646 0.013926 -0.075136 -0.001957 2 0.001077 0.001819 -0.006975 0.029570 0.000824 3 -0.003302 0.006161 0.000583 -0.023586 0.005854 4 0.006974 -0.015657 0.000219 0.016568 0.004597 5 0.005010 -0.008476 0.007476 0.047896 -0.005622 ... ... ... ... ... ... 2154 0.005785 0.004940 -0.023467 -0.014291 0.001349 2155 0.014971 0.036133 -0.001842 -0.007605 -0.011685 2156 0.006282 0.064420 0.004112 0.015402 0.004403 2157 0.002626 0.034169 0.003575 0.006245 -0.008100 2158 0.015374 -0.025590 0.001625 -0.015053 -0.001168 [2158 rows x 5 columns] We are now ready to run the correlation study. Let’s now subtract everything by the risk-free rate (dropping the rfr itself): risk_free_excess = returns.drop(columns="TBill").apply(lambda x: x-returns.TBill) risk_free_excess NYSE LMT TWTR MCD 1 -0.038846 -0.035720 -0.124783 -0.051603 2 -0.000742 -0.008794 0.027751 -0.000995 3 -0.009463 -0.005577 -0.029747 -0.000307 4 0.022630 0.015875 0.032225 0.020254 5 0.013486 0.015952 0.056372 0.002854 ... ... ... ... ... 2154 0.000845 -0.028406 -0.019231 -0.003591 2155 -0.021162 -0.037975 -0.043738 -0.047818 2156 -0.058138 -0.060308 -0.049017 -0.060017 2157 -0.031543 -0.030593 -0.027924 -0.042269 2158 0.040964 0.027215 0.010537 0.024422 [2158 rows x 4 columns] Actual Regression It is now time to perform the actual linear regression! import statsmodels.api as sm Let’s work with Lockheed Martin first, fitting an ordinary least squares. Remember that the OLS functions reads the endogenous variable first (for us, the return of the asset.) # add a column of ones to our input market excess returns nyse_with_bias = sm.add_constant(risk_free_excess.NYSE) # perform linreg lmt_model = sm.OLS(risk_free_excess.LMT, nyse_with_bias).fit() lmt_model.summary() OLS Regression Results ============================================================================== Dep. Variable: LMT R-squared: 0.859 Model: OLS Adj. R-squared: 0.859 Method: Least Squares F-statistic: 1.312e+04 Date: Mon, 31 Oct 2022 Prob (F-statistic): 0.00 Time: 10:39:24 Log-Likelihood: 6318.9 No. Observations: 2158 AIC: -1.263e+04 Df Residuals: 2156 BIC: -1.262e+04 Df Model: 1 Covariance Type: nonrobust ============================================================================== coef std err t P>|t| [0.025 0.975] ------------------------------------------------------------------------------ const 0.0004 0.000 1.311 0.190 -0.000 0.001 NYSE 0.9449 0.008 114.552 0.000 0.929 0.961 ============================================================================== Omnibus: 423.969 Durbin-Watson: 1.965 Prob(Omnibus): 0.000 Jarque-Bera (JB): 11575.074 Skew: -0.160 Prob(JB): 0.00 Kurtosis: 14.341 Cond. No. 29.6 ============================================================================== Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified. Based on the constants row, we can see that—within \(95\%\) confidence—the intercept is generally \(0\) and CAPM applies. However, we do see a slight positive compared to the market. Furthermore, we can see that the regression has a beta value of \(0.9449\) — according the CAPM model, it being slightly undervarying that the market. We can continue with the other stocks. # perform linreg mcd_model = sm.OLS(risk_free_excess.MCD, nyse_with_bias).fit() mcd_model.summary() OLS Regression Results ============================================================================== Dep. Variable: MCD R-squared: 0.887 Model: OLS Adj. R-squared: 0.887 Method: Least Squares F-statistic: 1.697e+04 Date: Mon, 31 Oct 2022 Prob (F-statistic): 0.00 Time: 10:39:24 Log-Likelihood: 6551.1 No. Observations: 2158 AIC: -1.310e+04 Df Residuals: 2156 BIC: -1.309e+04 Df Model: 1 Covariance Type: nonrobust ============================================================================== coef std err t P>|t| [0.025 0.975] ------------------------------------------------------------------------------ const 0.0003 0.000 1.004 0.315 -0.000 0.001 NYSE 0.9651 0.007 130.287 0.000 0.951 0.980 ============================================================================== Omnibus: 323.911 Durbin-Watson: 1.988 Prob(Omnibus): 0.000 Jarque-Bera (JB): 3032.550 Skew: 0.395 Prob(JB): 0.00 Kurtosis: 8.753 Cond. No. 29.6 ============================================================================== Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified. Same thing as before, we are within \(95\%\) confidence having a intercept of \(0\) (with a slight positive edge), and it looks like MacDonald’s vary a little bit more than Lockheed Martin. The food industry is probably a tougher business than that in defense. Lastly, to analyze the recently delisted Twitter! # perform linreg twtr_model = sm.OLS(risk_free_excess.TWTR, nyse_with_bias).fit() twtr_model.summary() OLS Regression Results ============================================================================== Dep. Variable: TWTR R-squared: 0.522 Model: OLS Adj. R-squared: 0.522 Method: Least Squares F-statistic: 2357. Date: Mon, 31 Oct 2022 Prob (F-statistic): 0.00 Time: 10:39:24 Log-Likelihood: 4307.1 No. Observations: 2158 AIC: -8610. Df Residuals: 2156 BIC: -8599. Df Model: 1 Covariance Type: nonrobust ============================================================================== coef std err t P>|t| [0.025 0.975] ------------------------------------------------------------------------------ const -0.0002 0.001 -0.346 0.730 -0.002 0.001 NYSE 1.0173 0.021 48.549 0.000 0.976 1.058 ============================================================================== Omnibus: 661.205 Durbin-Watson: 1.986 Prob(Omnibus): 0.000 Jarque-Bera (JB): 15925.609 Skew: -0.883 Prob(JB): 0.00 Kurtosis: 16.191 Cond. No. 29.6 ============================================================================== Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified. Evidently, Twitter is much more variable. It looks like it has a nontrivial bias (the intercept being -0.001 being within the \(95\%\) confidence band — that the security is possibly significantly underperforming the CAPM expectation in the market.) Furthermore, we have a positive beta value: that the asset is more variable than the market. manual regression We can also use the betas formula to manually calculate what we expect for the beta values (i.e. as if they were one IID random variable.) risk_free_cov = risk_free_excess.cov() risk_free_cov NYSE LMT TWTR MCD NYSE 0.001143 0.001080 0.001163 0.001103 LMT 0.001080 0.001188 0.001116 0.001083 TWTR 0.001163 0.001116 0.002264 0.001155 MCD 0.001103 0.001083 0.001155 0.001200 Finally, to construct the beta values. Recall that: \begin{equation} \beta_{im} := \frac{Cov[(R_{i}-R_{f}),(R_{m}-R_{f})]}{Var[R_{m}-R_{f}]} \end{equation} and that: \begin{equation} Var[X] = Cov[X,X], \forall X \end{equation} # get the market variance (covariance with itself) market_variation = risk_free_cov.NYSE.NYSE # calculate betas betas = {"LMT": (risk_free_cov.LMT.NYSE/market_variation), "TWTR": (risk_free_cov.TWTR.NYSE/market_variation), "MCD": (risk_free_cov.MCD.NYSE/market_variation)} # and make dataframe betas = pd.Series(betas) betas LMT 0.944899 TWTR 1.017294 MCD 0.965081 dtype: float64 Apparently, all of our assets swing less than the overall NYSE market! Especially Lockheed—it is only \(94.4\%\) of the market variation. Furthermore, it is interesting to see that Twitter swings much more dramatically compared to the market. Fund creation We will now create two funds with the three securities, one with equal parts and one which attempts to maximizes the bias (max returns) while minimizing the beta variance value compared to the market. “Equal-Parts Fund” (“Fund 1”) We will now create a fund in equal parts. Here it is: fund_1_returns = returns.LMT + returns.TWTR + returns.MCD fund_1_returns 1 -0.063167 2 0.023420 3 -0.017149 4 0.021384 5 0.049750 ... 2154 -0.036409 2155 -0.021132 2156 0.023917 2157 0.001720 2158 -0.014596 Length: 2158, dtype: float64 We will calculate the excess returns of this fund: fund_1_excess = fund_1_returns-returns.TBill fund_1_excess 1 -0.112813 2 0.021600 3 -0.023310 4 0.037041 5 0.058226 ... 2154 -0.041349 2155 -0.057265 2156 -0.040503 2157 -0.032449 2158 0.010994 Length: 2158, dtype: float64 And then perform a regression # perform linreg fund_1_model = sm.OLS(fund_1_excess, nyse_with_bias).fit() fund_1_model.summary() OLS Regression Results ============================================================================== Dep. Variable: y R-squared: 0.473 Model: OLS Adj. R-squared: 0.473 Method: Least Squares F-statistic: 1935. Date: Mon, 31 Oct 2022 Prob (F-statistic): 3.01e-302 Time: 10:39:24 Log-Likelihood: 3869.5 No. Observations: 2158 AIC: -7735. Df Residuals: 2156 BIC: -7724. Df Model: 1 Covariance Type: nonrobust ============================================================================== coef std err t P>|t| [0.025 0.975] ------------------------------------------------------------------------------ const 0.0007 0.001 0.841 0.401 -0.001 0.002 NYSE 1.1290 0.026 43.993 0.000 1.079 1.179 ============================================================================== Omnibus: 600.456 Durbin-Watson: 2.022 Prob(Omnibus): 0.000 Jarque-Bera (JB): 8416.514 Skew: -0.914 Prob(JB): 0.00 Kurtosis: 12.501 Cond. No. 29.6 ============================================================================== Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified. Surprisingly, we have now created a significantly riskier investment that, though riskier, generates a much higher probability of reward (\(+0.001\) is now within the \(99\%\) band!) A Better Fund To me, this is the payoff of this assignment. We will now use CAPM to create the “best” fund combination—given some variance, the funds which match CAPM. To do this, let’s create a generic linear combination of the assets. import sympy as sym x = sym.Symbol('x') y = sym.Symbol('y') z = sym.Symbol('z') fund_2_returns = x*returns.LMT + y*returns.TWTR + z*returns.MCD fund_2_returns 1 0.0139260753744255*x - 0.0751364261353569*y - ... 2 -0.00697525170622448*x + 0.0295704573211193*y ... 3 0.000583132897928884*x - 0.0235859990058791*y ... 4 0.000218587198947517*x + 0.016568426347233*y +... 5 0.00747599199607762*x + 0.0478955096700351*y -... ... 2154 -0.0234665578621085*x - 0.0142913301107561*y +... 2155 -0.00184214468578059*x - 0.0076045993852194*y ... 2156 0.00411172646842317*x + 0.0154024001854269*y +... 2157 0.00357547337231878*x + 0.0062445563228315*y -... 2158 0.00162509910496933*x - 0.0150529686289622*y -... Length: 2158, dtype: object Excellent. We will also calculate the excess returns of this fund: fund_2_excess = fund_2_returns-returns.TBill Y = fund_2_excess.to_numpy() Y [0.0139260753744255*x - 0.0751364261353569*y - 0.00195664549320096*z - 0.0496463208073039 -0.00697525170622448*x + 0.0295704573211193*y + 0.000824317408861575*z - 0.00181917459665826 0.000583132897928884*x - 0.0235859990058791*y + 0.00585367525146019*z - 0.00616055581298536 ... 0.00411172646842317*x + 0.0154024001854269*y + 0.00440300114913317*z - 0.0644196927849867 0.00357547337231878*x + 0.0062445563228315*y - 0.0081004573348249*z - 0.0341688956152497 0.00162509910496933*x - 0.0150529686289622*y - 0.00116834209450634*z + 0.0255902303732043] We cast this type to a numpy array because we are about to perform some matrix operations upon it. Now, let us perform the actual linear regression ourselves. Recall that the pseudoinverse linear regression estimator is: \begin{equation} \beta = (X^{T}X)^{-1}X^{T}Y \end{equation} We have already our \(Y\) as a vector (above), and our \(X\) is: X = nyse_with_bias.to_numpy() X [[ 1.00000000e+00 -3.88457302e-02] [ 1.00000000e+00 -7.42217926e-04] [ 1.00000000e+00 -9.46284244e-03] ... [ 1.00000000e+00 -5.81378271e-02] [ 1.00000000e+00 -3.15429207e-02] [ 1.00000000e+00 4.09643405e-02]] We now have our matrices, let’s perform the linear regression! linear_model = np.linalg.inv((X.transpose()@X))@X.transpose()@Y linear_model [0.000544056413840724*x - 6.62061061591867e-5*y + 0.000429966553373172*z - 0.000178620725465344 0.0457830563134785*x + 0.118178191274045*y + 0.0659651260604729*z + 0.899115719100281] Excellent. So we now have two rows; the top row represents the “bias”—how much deviation there is from CAPM, and the bottom row represents the “rate”—the “beta” value which represents how much excess variance there is. We can will solve for a combination of solutions to give us specific values of returns vs risk. We will set the asset to learn exactly as much as the market (i.e. no bias). deviance_expr = linear_model[0] deviance_expr 0.000544056413840724*x - 6.62061061591867e-5*y + 0.000429966553373172*z - 0.000178620725465344 We will now try to make variance exactly as much as that in the market. risk_expr = linear_model[1] - 1 risk_expr 0.0457830563134785*x + 0.118178191274045*y + 0.0659651260604729*z - 0.100884280899719 Let us now calculate the boundary condition of our optimization problem by solving an expression in these two expressions. solution = sym.solvers.solve([deviance_expr, risk_expr], x,y,z) solution {x: 0.412737013327711 - 0.819584899551304*z, y: 0.693765220909132 - 0.24067066980814*z} Excellent. Let us recalculate our optimization objective (“deviance”—return) in terms of these new solutions. We aim now to maximize this expression by minimizing (i.e. our optimizer minimizes) the negative thereof—recalling that scypy works as a minimizer. optim_objective = deviance_expr.subs(solution)-1e2 optim_objective -5.04831636563563e-19*z - 100.0 We can now use this value to solve for a \(z\) value. optim_solution = sym.solvers.solve([optim_objective], z) optim_solution {z: -1.98085842402250e+20} Excellent. We can now solve for the rest of our values. z0 = float(optim_solution[z]) x0 = solution[x].subs(z, z0) y0 = solution[y].subs(z, z0) (x0,y0,z0) (1.62348165247784e+16, 4.76734523704593e+15, -1.980858424022502e+16) This would create the following plan: # solution fund_2_nobias = x0*returns.LMT + y0*returns.TWTR + z0*returns.MCD fund_2_nobias.mean() 0.009168283711770158 Recall that this is the performance of the balanced portfolio: fund_1_returns.mean() 0.0009224705380695683 Finally, let’s plot the prices of our various funds: import matplotlib.pyplot as plt import matplotlib.dates as mdates import seaborn as sns from datetime import datetime sns.set() fund_2_price = x0*df.LMT + y0*df.TWTR + z0*df.MCD fund_1_price = df.LMT + df.TWTR fund_l_price = df.LMT fund_t_price = df.TWTR dates = df.Date.apply(lambda x:datetime.strptime(x, "%m/%d/%Y %H:%M:%S")) sns.lineplot(x=dates, y=fund_2_price.apply(sym.Float).astype(float)) sns.lineplot(x=dates, y=fund_1_price.apply(sym.Float).astype(float)) sns.lineplot(x=dates, y=fund_l_price.apply(sym.Float).astype(float)) sns.lineplot(x=dates, y=fund_t_price.apply(sym.Float).astype(float)) plt.gca().xaxis.set_major_locator(mdates.YearLocator()) plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y')) plt.gca().set_ylabel("Price") plt.show() None Recall that we didn’t actually buy any MacDonald’s. So, we have—blue, being our “optimal” portfolio, yellow, our balanced portfolio, green, being Lockheed, and red, being Twitter. Our portfolio works surprisingly well!