[This article was first published on R – insightR, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

By Gabriel Vasconcelos

## Overview

There are several ways to do portfolio optimization out there, each with its advantages and disadvantages. We already discussed some techniques here. Today I am going to show another method to perform portfolio optimization that works very well in large datasets because it produces very robust weights, which results in a good out-of-sample performance. This technique is called Parametric Portfolio Policies (PPP) and it was proposed by Brandt, Santa-Clara and Valkanov in 2009 (click here to read the full article).

A portfolio is a set of weights $w_{i,t}$ that states the proportion of wealth in asset $i$ as time $t$. For example, in a Markowitz portfolio we get $w_{i,t}$ by minimizing the portfolio variance subject to some target return. This is like estimating one weight for each asset if you think of weights as parameters of a model. If you have 1000 assets you will have to estimate 1000 parameters. Naturally, you can’t expect a good out-of-sample performance from a simple Markowitz portfolio with 1000 assets just as you can’t expect a good performance from a simple regression with 1000 variables. The Parametric Portfolio Policy (PPP) is made exactly to this high-dimensional case. It parametrizes the weights on some asset characteristics in a way that we have to estimate only one parameter for each characteristic. Consider the following portfolio policy:

$\displaystyle w_{i,t} = \bar{w}_{i,t} + \frac{1}{N_t}\theta'\hat{x}_{i,t}$

Here, the portfolio $w_{i,t}$ is a function of a base portfolio $\bar{w}_{i,t}$, which will be a simple equal weighted portfolio in my example. $\hat{x}_{i,t}$ is a vector of characteristics for asset $i$ at time $t$, $\theta$ is the vector of parameters and $N_t$ is the number of available assets. Note that $\hat{x}_{i,t}$ is standardized cross-sectionally to make sure the sum of the weights in each period is one. If we have two characteristics we only need to estimate two parameters and use the policy equation to recover all weights in the portfolio. Unfortunately using a policy like this makes the portfolio much more restricted than the Markowitz case. In a Markowitz portfolio weights can be anything as long as the sum 1 but in the PPP the portfolio must satisfy the policy equation.

Now that we have the policy we need a way to obtain $\theta$. We can use utility functions and try to maximize the portfolio utility. In the equation below we have the populational case in the left and the sample case in the right.

$\displaystyle \max_{\theta}E[u(r_{p,t+1})] \sim \max_{\theta} \frac{1}{T}\sum_{t=0}^{T-1}u(r_{p,t+1})$

where:

$\displaystyle r_{p,t+1} = \sum_{i=1}^{N_t} w_{i,t}r_{i,t+1}$

Putting all together we have:

$\displaystyle \max_{\theta} \frac{1}{T}\sum_{t=0}^{T-1}u\left(\left( \bar{w}_{i,t} + \frac{1}{N_t}\theta'\hat{x}_{i,t}\right)r_{i,t+1}\right)$

where the utility function will be the Constant Relative Risk Aversion (CRRA) with relative risk aversion $\gamma$, define as:

$\displaystyle u(x)=\frac{x^{1-\gamma}-1}{1-\gamma}$

## The data

The PPP parametrizes the weights on characteristics of the stocks, and although we can get a lot of characteristics from returns such as momentum and variance, we are likely to need something else to have good results. In the original paper Brandt Santa-Clara and Valkanov used momentum, book-to-market ratio and market capitalization (mktcap) . I will use only a momentum measure and the mktcap because they are easier to get and need no complex treatment.

The code below uses the package tidyquant to download data from all S&P500 stocks from Yahoo Finance. The code takes some time to download all data. If you want to replicate the exact example I will show you can also download the data from here. My data starts at 2008-01-02 and ends at 2018-02-06. If you use a different time you will obviously find different results.

library(tidyverse)
library(reshape2)
library(lubridate)
library(tidyquant)

prices=tq_index("sp500") %>%
tq_get(get = "stock.prices", complete_cases = TRUE)

## # A tibble: 6 x 12
##   symbol company  weight sector  share… date        open  high   low close
##   <chr>  <chr>     <dbl> <chr>    <dbl> <date>     <dbl> <dbl> <dbl> <dbl>
## 1 MSFT   Microso… 0.0306 Inform… 9.22e⁷ 2008-01-02  45.9  46.1  44.9  35.2
## 2 MSFT   Microso… 0.0306 Inform… 9.22e⁷ 2008-01-03  45.2  45.7  44.7  35.4
## 3 MSFT   Microso… 0.0306 Inform… 9.22e⁷ 2008-01-04  45.1  45.2  43.7  34.4
## 4 MSFT   Microso… 0.0306 Inform… 9.22e⁷ 2008-01-07  44.3  44.6  43.9  34.6
## 5 MSFT   Microso… 0.0306 Inform… 9.22e⁷ 2008-01-08  44.5  44.5  42.8  33.5
## 6 MSFT   Microso… 0.0306 Inform… 9.22e⁷ 2008-01-09  42.8  44.3  42.8  34.4
## # ... with 2 more variables: volume <dbl>, adjusted <dbl>


I will build the portfolio optimization on the adjusted prices. We also need to select close prices and shares_held to calculate the market capitalization. Additionally, I will optimize the portfolio on monthly returns and the data must be aggregated.

## == Treat Data == ##
prices = prices %>%
select(symbol, date, adjusted, close ,shares_held) %>%
group_by(year(date), month(date), symbol) %>%
summarize( adjusted = tail(adjusted, 1),
close = tail(close, 1),
mktcap = tail(close * shares_held, 1)) %>%
mutate(date = make_date(year(date), month(date)))


Next we need to transform the data into time-series panel with the function acast and calculate log-returns from the adjusted prices. Note that the characteristics in $t$ are used to build a portfolio for $t+1$. Therefore, the characteristics must be lagged one month from the returns. I also removed stocks with NAs. The PPP works fine in unbalanced panels but I will use them balanced to keep it simple. Just keep in mind that if you use it in unbalanced panels $N_t$ and the base portfolio $\bar{w}_{i,t}$ must be adjusted.

## == Reshape variables to matrices == ##
# = dropped first observatio in the mktcap because we = #
# = loose one observatio to calculate the returns     = #
mktcap=acast(prices,date~symbol,value.var = "mktcap")[-1,]
prices=acast(prices,date~symbol,value.var = "adjusted")

## == Calculate log returns == ##
ret=diff(log(prices))

## == Remove Stocks with NAs == ##
drop = which( is.na(colSums(ret)) | is.na(colSums(mktcap)) )
ret = ret[, - drop]
mktcap = mktcap[ , - drop]

## == Remove last observation of mktcap because it has == ##
## == to be lagged one month from the returns          == ##
mktcap = mktcap[- nrow(mktcap), ]


The next step is to generate a variable that is the cumulative return of the last 12 months that will be used as a momentum characteristic. We loose some observations in this procedure and the rest of the data must be adjusted. The characteristics must also be standardized cross-sectionally to have mean zero and variance one. The next chunk also defines the risk aversion coefficient as 5. Feel free to play with the risk aversion. If you increase it to much the portfolio will be the same as the equal weighted and if you decrease the risk aversion you are likely to have portfolios with more variance and leverage.

## == Generate 12 month cumulative returns  == ##
## == I used sub because returns are in log == ##
m12 = matrix(NA, nrow(ret), ncol(ret))
for(i in 13:nrow(m12)){
m12[i, ] = apply(ret[(i - 12):(i - 1), ],2 , sum )
}

## == Remove last 12 observations because we == ##
## == don't the first 12 m12.                == ##
ret = ret[-c(1:12), ]
m12 = m12[-c(1:12), ]
# = only 11 in the mktcap because we already = #
# = removed one                              = #
mktcap = mktcap[-c(1:11), ]

## == This is the weigth of each stock in a == #
## == equal weighted portfolio              == #
nt = wb = 1/ncol(ret)

## = Risk aversion == ##
rr = 5

## == Some plots == ##
par(mfrow = c(2, 2))
matplot(ret[ ,1:10], type = "l", ylab = "monthly returns")
matplot(m12[ ,1:10], type = "l", ylab = "12 month returns")
matplot(mktcap[ ,1:10], type = "l", ylab = "market captalization")

## == Characteristics must be standardized == ##
## == on the crossection == ##
m12=t(apply(m12,1,scale))
mktcap=t(apply(mktcap,1,scale))


## Optimizing

Now let’s define the optimization function and use it to find our portfolio. I will use the function in 60 expanding windows to simulate a portfolio that is re-estimated and updated every month. The interpretation is that in the end of every month we optimize with all the data we have until then and build our portfolio for the next month.

## == Portfolio Optimization Function == ##
pps = function(x, wb, nt, ret, m12, mktcap, rr){
wi = wb + nt * (x[1]*m12 + x[2]*mktcap)
wret = rowSums(wi*ret)
ut = ((1 + wret)^(1 - rr))/(1 - rr)
u = -mean(ut)
return(u)
}

## == Loop the function in 60 expanding windows == ##
res_save = matrix(NA, 60, 2)
weights = matrix(NA, 60, ncol(ret))
for(i in 1:60){
opt = optim(c(0, 0), pps, wb = wb, nt = nt, ret = ret[1:(48 + i), ],
m12 = m12[1:(48 + i), ], mktcap = mktcap[1:(48 + i), ],
rr = rr, method = "BFGS")
res = opt$par w = wb + nt*(res[1]*m12[i + 49, ] + res[2]*mktcap[i + 49, ]) weights[i, ] = w res_save[i, ] = res }  The code below calculates the out-of-sample cumulative returns of our optimized portfolio and the equal weighted (EW) portfolio to use as benchmark. I also calculated the returns of a look-ahead portfolio with the best 100 stocks in the S&P500 during our period of analysis. The code ends with a plot comparing the cumulative returns of the tree portfolios. ## == Calculate out-of-sample returns == ## ret_fit = rowSums(weights*tail(ret, 60)) ret_EW = rowSums(nt*tail(ret, 60)) ## == Accumulate out-of-sample returns == ## acc_fit = cumsum(ret_fit) acc_EW = cumsum(ret_EW) ## == Calculate top 100 portfolio == ## acc_ret = apply(ret[-c(1:49), ], 2, cumsum) top100 = order(acc_ret[nrow(acc_ret), ] ,decreasing = TRUE)[1:100] acc_top100 = rowMeans(acc_ret[ ,top100]) df = data.frame(date = tail(as.Date(rownames(ret)), 60), opt = acc_fit,EW = acc_EW,top100 = acc_top100) dfm = melt(df, id.vars = "date") colnames(dfm) = c("date", "portfolio", "return") ggplot(dfm) + geom_line(aes(x = date, y = return, color = portfolio))  As you can see, our optimal portfolio is far superior than the EW portfolio and it beats the top100 portfolio most of the time. The next plot shows the parameters of the PPP in the 60 expanding windows. In general, our policy is to buy more of winners (high momentum) and buy less of bigger companies. df = data.frame(date = tail(as.Date(rownames(ret)), 60), res_save) colnames(df)=c("date","m12","mktcap") dfm = melt(df, id.vars = "date") colnames(dfm) = c("date", "characteristic", "parameter") ggplot(dfm) + geom_line(aes(x = date, y = parameter, color = characteristic)) + facet_wrap(~ characteristic, scales = "free_y", ncol = 1)  ## Important Considerations The question now is: does the optimal portfolio I created generate realistic results??? The answer is NO! The results were very impressive, the PPP was able to beat a look-ahead portfolio, which is very unlikely to happen. The cause of this super performance is the lack of restrictions on short positions. Our portfolio assumes that we have infinite credit at no cost to borrow money and invest in stocks. Look at the plot below: lv = apply(weights, 1, function(x) sum(x[x<0])) plot(abs(lv),type="l")  This plot shows the leverage we had in each month in the out-of-sample period, and it is a lot. In some months we had 4.5 times our wealth of borrowed money. Fortunately, it is very easy to put some restrictions in the PPP. All we have to do is add a penalty in the objective function that penalizes short positions when they are bigger than some target we choose. The function below penalizes leverage values bigger than two. The code chunk ends with the leverage comparison plot, which shows that the penalization kept the leverage very close to our target of two. ## == Model with leverage penalty == ## pps1 = function(x, wb, nt, ret, m12, mktcap, rr){ wi = wb + nt * (x[1]*m12 + x[2]*mktcap) wret = rowSums(wi*ret) ut = ((1 + wret)^(1 - rr))/(1 - rr) u = -mean(ut) lv = abs(apply(wi, 1, function(x) sum(x[x<0]))) pen=abs(sum((lv>2)*lv)) return(u+pen) } ## == Loop the function in 60 expanding windows == ## res_save1 = matrix(NA, 60, 2) weights1 = matrix(NA, 60, ncol(ret)) for(i in 1:60){ opt = optim(c(0, 0), pps1, wb = wb, nt = nt, ret = ret[1:(48 + i), ], m12 = m12[1:(48 + i), ], mktcap = mktcap[1:(48 + i), ], rr = rr, method = "BFGS") res = opt$par
w = wb + nt*(res[1]*m12[i + 49, ] + res[2]*mktcap[i + 49, ])
weights1[i, ] = w
res_save1[i, ] = res
}

## == Calculate out-of-sample returns == ##
ret_fit1 = rowSums(weights1*tail(ret, 60))

## == Accumulate out-of-sample returns == ##
acc_fit1 = cumsum(ret_fit1)

## == Leverage analysis == ##
lv1 = apply(weights1, 1, function(x) sum(x[x<0]))
df=data.frame(date = tail(as.Date(rownames(ret)), 60),lv = abs(lv), lv_pen = abs(lv1))
dfm = melt(df, id.vars = "date")
colnames(dfm) = c("date", "model", "leverage")
ggplot(dfm) + geom_line(aes(x = date, y = leverage, color = model))


However, the inclusion of short selling restrictions comes at a cost. Look at the plot below! The returns in the penalizes portfolio are smaller, which is natural because we now have a more restricted model.

df = data.frame(date = tail(as.Date(rownames(ret)), 60),
opt = acc_fit, opt_pen=acc_fit1, EW = acc_EW,top100 = acc_top100)
dfm = melt(df, id.vars = "date")
colnames(dfm) = c("date", "portfolio", "return")
ggplot(dfm) + geom_line(aes(x = date, y = return, color = portfolio))


## References

Brandt, Michael W., Pedro Santa-Clara, and Rossen Valkanov. “Parametric portfolio policies: Exploiting characteristics in the cross-section of equity returns.” The Review of Financial Studies 22.9 (2009): 3411-3447.

Medeiros, Marcelo, Artur M Passos, and Gabriel FR Vasconcelos. “Parametric Portfolio Selection: Evaluating and Comparing to Markowitz Portfolios.” Revista Brasileira de Finanças 12.2 (2014).

To leave a comment for the author, please follow the link and comment on their blog: R – insightR.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

# Never miss an update! Subscribe to R-bloggers to receive e-mails with the latest R posts.(You will not see this message again.)

Click here to close (This popup will not appear again)