Computing market price of risk for more than one risky asset

[Clem]

Hi everybody,

I've written a matlab code in order to calculate the market price of risk for historic data. But the market price of risk with more than one asset is somehow not realistic.

Just to have a commen understanding:

market price of risk = (drift - risk-free rate)/volatility

Now, for one asset this is pretty simple. I calculate the mean annual drift rate and the annual volatility (by using log-returns). With an average risk-free rate of about 4% I get a market price of risk of about 0.40. Of course the result depends on the assset and the time horizon of my financial data but the result is usually between 0.3 and 0.5. Seems plausible to me (hope you agree).

However, when I calculate the market price of risk with several assets I run into some troubles.

My calculations:

1. calculate the log-returns of each asset over time
2. calculate variance-covariance matrix (vcm) by using the log-returns
3. calculate volatility matrix by performing the cholesky decomposition (chol(vcm))
4. annualize the volatility matrix
5. calculate mean annual drift-rates for each asset

up to that point the data looks quite plausible. I have volatilities between 10 and 25%, covariances from -2% to +5% and drift rates of 7% to 12%.

Now I would like to calculate the market price of risk of all the assets as a vector. Of course I have to take the covariances into account, so I can not simply calculate it for each asset isolated.

Assuming that I have n assets, I calculated it as follows:
market price of risk vector = (drift - risk-free)/volatility,
where drift is a nx1 vector with the mean drift-rates, risk-free is a nx1 vector with the risk-free rate and volatility is the volatility matrix (calculated in step 4, see above).

The troubling thing is the result. I get a nx1 vector with some components below -1 and above +1. This seems to be not really plausible.

Just a side note: I can not calculate the market price of risk as an aggregated value (e.g. by aggregating the assets with equal weights into an index) because I need the vector for portfolio optimization later on.

Has anyone done a similar calculations and can give me some advice? Every comment is appreciated.
Thank you.

 

[TheComplexUnit]

Greetings,

Whilst I have not performed an empirical analysis of this problem, I suppose I can add theoretical support to your algorithm. I say suppose, considering that there are quite a few details which (understandably) are glossed over.

In order to see if our approaches indeed are in agreement, I here offer my own description of the solution approach:

Suppose we have n different traded financial securities \(\{S_1, S_2,...,S_n \}\) which form a complete market. For each asset, \(S_{i}\), we have a time series of m+1 consecutive observations: \(\{S_{i,t}\}_{t=1}^{m+1}\). Finally, we imagine that each asset follows (scalar) Geometric Brownian Motion \(dS_i = S_i [\mu_i dt + \sigma_i dW_i ]\) where \(dW_i dW_j = \rho_{i,j} dt \).

In order to estimate the various quantities "drift", "diffusion" and "correlation" we use standard statistical results. Specifically, recalling that \(R_{i,t} \equiv \ln (S_{i,t}/S_{i,t-1}) \), where \(t = 2,3,...,m+1\), is distributed normally with a mean of \((\mu_i-\tfrac{1}{2} \sigma_i^2)\Delta t\) and a variance of \(\sigma_i^2 \Delta t\) we may realise that the covariance matrix is estimated as

\(\boldsymbol{\Sigma}_{i,j} = \text{cov}[R_i,R_j] = \rho_{i,j} \sigma_i \sigma_j \Delta t \approx \frac{1}{m-1} \sum_{t=2}^{m+1} (R_{i,t} - \bar{R}_i)(R_{j,t}-\bar{R}_j)\).

where \(\bar{R}_i \approx \frac{1}{m} \sum_{t=2}^{m+1} R_{i} \). Backing out \(\mu_i, \sigma_i\) and \(\rho_{i,j}\) \(\forall i \forall j\) from these considerations is (here) deemed self-evident (as you clearly suggest, the importance lies in keeping track of the \(\Delta t\)).

Now to calculate the market price of risk vector we must first rewrite the governing dynamics in terms of a standard Wiener process, i.e. a Wiener vector where all components are independent. Presently, what we have is the vector equation \(d \boldsymbol{S} = \text{diag}[\boldsymbol{S}][\boldsymbol{\mu}dt + \text{diag}[\boldsymbol{\sigma}] d\boldsymbol{W}]\) where \(\boldsymbol{W} \equiv (W_1,W_2,...,W_n)\) consists of correlated components. Again, as you clearly indicate, the answer lies in the Cholesky decomposition. Specifically, since \(\text{diag}[\boldsymbol{\sigma}] d\boldsymbol{W} \sim N(\boldsymbol{0}, \boldsymbol{\Sigma}) \) then \(\text{diag}[\boldsymbol{\sigma}] d\boldsymbol{W}\) may be rewritten as \(\boldsymbol{L} \boldsymbol{Z}\), where \(\boldsymbol{L}\) is the lower triangular matrix in \(\boldsymbol{\Sigma} = \boldsymbol{L} \boldsymbol{L}^T\) and \(\boldsymbol{Z} \sim N(\boldsymbol{0},\boldsymbol{I})\). In more standard form, this means that our dynamics may be rewritten as \(d \boldsymbol{S} = \text{diag}(\boldsymbol{S})[\boldsymbol{\mu}dt + \bar{\boldsymbol{\sigma}} d\bar{\boldsymbol{W}}]\) where \(\bar{\boldsymbol{\sigma}} \equiv \boldsymbol{L} (dt)^{-0.5}\) and \(\bar{\boldsymbol{W}}\) is a standard Wiener process.

The market price of risk vector, \(\boldsymbol{\lambda}\) may now be expressed (per definitionem) as

\(\boldsymbol{\lambda} \equiv \bar{\boldsymbol{\sigma}}^{-1} (\boldsymbol{\mu}-r \boldsymbol{1}) \)

where r is the risk free rate.

 

[Clem]

Thanks for your input. I think your theoretic elaboration is inline with my approach. I calculated drift and diffusion by using log returns. The calculation of the volatility matrix is then straight forward since matlab has a build in function for chlolesky decomposition.

 

[themarketprofessor]

Greetings,

Whilst I have not performed an empirical analysis of this problem, I suppose I can add theoretical support to your algorithm. I say suppose, considering that there are quite a few details which (understandably) are glossed over.

In order to see if our approaches indeed are in agreement, I here offer my own description of the solution approach:

Suppose we have n different traded financial securities \(\{S_1, S_2,...,S_n \}\) which form a complete market. For each asset, \(S_{i}\), we have a time series of m+1 consecutive observations: \(\{S_{i,t}\}_{t=1}^{m+1}\). Finally, we imagine that each asset follows (scalar) Geometric Brownian Motion \(dS_i = S_i [\mu_i dt + \sigma_i dW_i ]\) where \(dW_i dW_j = \rho_{i,j} dt \).

In order to estimate the various quantities "drift", "diffusion" and "correlation" we use standard statistical results. Specifically, recalling that \(R_{i,t} \equiv \ln (S_{i,t}/S_{i,t-1}) \), where \(t = 2,3,...,m+1\), is distributed normally with a mean of \((\mu_i-\tfrac{1}{2} \sigma_i^2)\Delta t\) and a variance of \(\sigma_i^2 \Delta t\) we may realise that the covariance matrix is estimated as

\(\boldsymbol{\Sigma}_{i,j} = \text{cov}[R_i,R_j] = \rho_{i,j} \sigma_i \sigma_j \Delta t \approx \frac{1}{m-1} \sum_{t=2}^{m+1} (R_{i,t} - \bar{R}_i)(R_{j,t}-\bar{R}_j)\).

where \(\bar{R}_i \approx \frac{1}{m} \sum_{t=2}^{m+1} R_{i} \). Backing out \(\mu_i, \sigma_i\) and \(\rho_{i,j}\) \(\forall i \forall j\) from these considerations is (here) deemed self-evident (as you clearly suggest, the importance lies in keeping track of the \(\Delta t\)).

Now to calculate the market price of risk vector we must first rewrite the governing dynamics in terms of a standard Wiener process, i.e. a Wiener vector where all components are independent. Presently, what we have is the vector equation \(d \boldsymbol{S} = \text{diag}[\boldsymbol{S}][\boldsymbol{\mu}dt + \text{diag}[\boldsymbol{\sigma}] d\boldsymbol{W}]\) where \(\boldsymbol{W} \equiv (W_1,W_2,...,W_n)\) consists of correlated components. Again, as you clearly indicate, the answer lies in the Cholesky decomposition. Specifically, since \(\text{diag}[\boldsymbol{\sigma}] d\boldsymbol{W} \sim N(\boldsymbol{0}, \boldsymbol{\Sigma}) \) then \(\text{diag}[\boldsymbol{\sigma}] d\boldsymbol{W}\) may be rewritten as \(\boldsymbol{L} \boldsymbol{Z}\), where \(\boldsymbol{L}\) is the lower triangular matrix in \(\boldsymbol{\Sigma} = \boldsymbol{L} \boldsymbol{L}^T\) and \(\boldsymbol{Z} \sim N(\boldsymbol{0},\boldsymbol{I})\). In more standard form, this means that our dynamics may be rewritten as \(d \boldsymbol{S} = \text{diag}(\boldsymbol{S})[\boldsymbol{\mu}dt + \bar{\boldsymbol{\sigma}} d\bar{\boldsymbol{W}}]\) where \(\bar{\boldsymbol{\sigma}} \equiv \boldsymbol{L} (dt)^{-0.5}\) and \(\bar{\boldsymbol{W}}\) is a standard Wiener process.

The market price of risk vector, \(\boldsymbol{\lambda}\) may now be expressed (per definitionem) as

\(\boldsymbol{\lambda} \equiv \bar{\boldsymbol{\sigma}}^{-1} (\boldsymbol{\mu}-r \boldsymbol{1}) \)

where r is the risk free rate.

One hard assumption is that all stocks are Lognormal processes with constant drift and vol. Hmmm...people can live with assumption? i guess so.

 

[Clem]

I'm just wondering, what is from a theoretic point of view the correct approach for calculating the volatility matrix: cholesky decomp. of the covariance matrix or the chlesky decomp of the correlation matrix?

thanks.

 

[TheComplexUnit]

I'm just wondering, what is from a theoretic point of view the correct approach for calculating the volatility matrix: cholesky decomp. of the covariance matrix or the chlesky decomp of the correlation matrix?
thanks.

You can do both.

Recall that in the dynamics \(d \boldsymbol{S} = \text{diag}[\boldsymbol{S}][\boldsymbol{\mu}dt + \text{diag}[\boldsymbol{\sigma}] d \boldsymbol{W}]\), the diffusion component \(\text{diag}[\boldsymbol{\sigma}] d \boldsymbol{W}\) has distribution \(N(\boldsymbol{0},\boldsymbol{\Sigma})\) where \(\boldsymbol{\Sigma}\) is the covariance matrix "of the log returns over the instant \([t,t+dt]\)" (i.e. with components \(\rho_{i,j} \sigma_i \sigma_j dt\)).

Alternatively, we could say that the diffusion component \(d \boldsymbol{W}\) has distribution \(N(\boldsymbol{0},\boldsymbol{\rho})\) where \(\boldsymbol{\rho}\) is the correlation matrix "of the log returns over the instant \([t,t+dt]\)" (i.e. with components \(\rho_{i,j} dt\)).

Now, our task is to rewrite the dynamics in terms of a standard Wiener process.

We can do this as follows: CD the covariance matrix \(\boldsymbol{\Sigma}\) as \(\boldsymbol{L}\boldsymbol{L}^T\) in which case \(\text{diag}[\boldsymbol{\sigma}] d \boldsymbol{W}\) is identically distributed to (can be replaced by) \(\boldsymbol{L}\boldsymbol{Z}\). Of course, \(\boldsymbol{Z}\) is not a standard Wiener process, but \(\sqrt{dt} \boldsymbol{Z}\) is. Therefore, \(\text{diag}[\boldsymbol{\sigma}] d \boldsymbol{W}\) is identically distributed to \((\boldsymbol{L}/\sqrt{dt}) d \bar{\boldsymbol{W}}\) where \(\bar{\boldsymbol{W}}\) is a standard Wiener process.

Alternatively, we could CD the correlation matrix \(\boldsymbol{\rho}\) as \(\boldsymbol{L}_{*}\boldsymbol{L}_{*}^T\) in which case \(d \boldsymbol{W}\) is identically distributed to (can be replaced by) \(\boldsymbol{L}_{*}\boldsymbol{Z}\). Again, \(\boldsymbol{Z}\) is not a standard Wiener process, but \(\sqrt{dt} \boldsymbol{Z}\) is. Therefore, \( d \boldsymbol{W}\) is identically distributed to \((\boldsymbol{L}_{*}/\sqrt{dt}) d \bar{\boldsymbol{W}}\) where \(\bar{\boldsymbol{W}}\) is a standard Wiener process.

Obviously, there's a relation between the two: specifically, since \(\boldsymbol{\Sigma}= \text{diag}[\boldsymbol{\sigma}] \boldsymbol{\rho} \text{diag}[\boldsymbol{\sigma}]\) then \(\boldsymbol{L} = \text{diag}[\boldsymbol{\sigma}] \boldsymbol{L}_{*}\).