Estimating Future Price Ranges with Volatility Cones

# R Calculate future price range and plot the result


# R Calculate future price range and plot the result












First I want to say that I've read this post (How to calculate future distribution of price using volatility?) but it doesn't help much.

Here is what I'm trying to do (values are not real)

Let's pretend that the last price for this security is 20.85 and what I'm trying to do is to determine, based on volatility, the range where the price should be in n days with X % of probability.

On my plot let's pretend that the blue area is a probability of 41.5% and the two pruple dots are on 20 days in the future. I'm looking for a function that would take a daily end of day price time serie, a probability, and a number of day and return a price range. For instance based on my plot:

- In 10 days the price should be between 19.65 and 20.05 with 93.4% chances

- In 20 days the price should be between 20.20 and 21.50 with 41.5% chances

I would like to be able to do that with a function like this :

```
future_price_range(eod_prices_time_serie, nb_of_days, probability)

#with the first exemple it would look like this:
future_price_range(pcef_eod, 10, 0.934)

# And return
[1] 19.65 20.05
```

Based on the post I'm refering to at the top of this question I already managed to do something but I'm not sure if I'm doing it right:

```
library(quantmod)
library(PerformanceAnalytics)

# Retrieving daily prices from yahoo finance 
pcef <- getSymbols("PCEF")

# Selecting only the column of adjusted prices
pcef_ad <- Ad(PCEF)
```

In the post I've linked at the top of my question it is said : "The distribution of the log of a stock price in n days is a normal distribution with mean of log(currentprice) ..."

```
# So i take the last price
tail(pcef_ad, 1)
# PCEF.Adjusted
# 2017-11-22         23.62

# And a convert it to log 
mean_last_price_log <- log(23.62)
```

" ... and standard deviation of volatility∗(√n/365.2425) if you're using calendar days, and assuming no dividends and 0% risk-free interest rate."

```
# I calculate daily returns 
pcef_daily_return <- dailyReturn(pcef_ad)

# To calculate the standard dev of returns
sd(pcef_daily_return)
# [,1]
# 0.005932502
```

Since I'm using the volatility of daily returns I guess I should not divide n by 365.2425. Am I right about this ?

```
# Let's say we want the price range in 10 days
n <- 10

sd_pcef_daily_return <- 0.005932502 * sqrt(n)

lg_dist <- rnorm(n = 10000, mean = mean_last_price_log, sd = sd)

# I convert it back to USD
us_dist <- exp(log_dist)
hist(us_dist)
```

I'm not sure what number should i take for the n parameter in rnorm. I took 10000 because it seems "enough" but what is the rule here ? Here is what i get:

```
### making a function of this

future_price_range <- function(eod_prices, days, probability){
  last_price <- tail(eod_prices,1)
  mean_log_price <- log(last_price)
  sd_daily_returns <- sd(dailyReturn(eod_prices)) * sqrt(days)
  log_dist <- rnorm(n = 10000, mean = mean_log_price, sd = sd_daily_returns)
  usd_dist <- exp(log_dist)
  return(hist(usd_dist))
}

future_price_range(eod_prices = pcef_ad, days = 100)
```

So with 100 days the probability of having a higher/lower price increases as expected:

But now I'm struggling to implement the "probability" parameter in my function and return a vector of the 2 prices for the range.

Could you confirm that I haven't done any mistake in my code and do you have any clues to help me to finish It.

It would also be great if i could output the same kind of plot I used at the beginning of this post. I guess I need to compute the price range for each days until the target date and use the data to plot it on the graph but I have no idea how to do that.

There is one last thing I'd like to do: I take a end of day daily price time series, i set a number of days and i get a price range for a set of probabilities:

```
function(price_data, nb_of_days)

99%    55.12    20.90
95%    54.34    21.36
90%    53.26    22.35
80%    49.78    24.12
...    ...      ...
10%    ...      ...
5%     ...      ...
1%     ...      ...
```

Any help would be much appreciated. Thank you so much in advance.

### EDIT : Answer from @vanguard2k (Thank you so much !)

I made a few adjustment from @vanguard2k answer to make the code reproductible. For more details, have a look at the post below

```
library(PerformanceAnalytics)
library(quantmod)
library(ggplot2)
library(lubridate)

# Retrieving price data from yahoo
pcef <- getSymbols("PCEF", auto.assign = FALSE)

# Selecting only adjusted EOD prices
pcef_ad <- Ad(pcef)

# tail(pcef_ad)
# PCEF.Adjusted
# 2017-11-15      23.21980
# 2017-11-16      23.40866
# 2017-11-17      23.51800
# 2017-11-20      23.53000
# 2017-11-21      23.57000
# 2017-11-22      23.62000

# Computing daily returns
pcef_daily_returns <- dailyReturn(pcef_ad)

# tail(pcef_daily_returns)
# daily.returns
# 2017-11-15 -0.0004278567
# 2017-11-16  0.0081334892
# 2017-11-17  0.0046709639
# 2017-11-20  0.0005102900
# 2017-11-21  0.0016999149
# 2017-11-22  0.0021213831

# Annualized Standard Deviation of daily returns
sigma <- sd(pcef_daily_returns)*sqrt(365)

# Why do we set a mean of zero ???
mean <- 0
probability <- 0.95

# Last observed price as "S"
S <- coredata(tail(pcef_ad,1))[1,1]

# Number of days to "forecast" as "n"
n <- 10

upper.bounds <- qnorm(probability,mean=mean*(0:n)/365,sd=sigma*sqrt((0:n)/365))
# [1] 0.000000000 0.009768138 0.013814233 0.016918911 0.019536276
# [6] 0.021842220 0.023926954 0.025844064 0.027628466 0.029304414
# [11] 0.030889564

lower.bounds <- qnorm(1-probability,mean=mean*(0:n)/365,sd=sigma*sqrt((0:n)/365))
# [1]  0.000000000 -0.009768138 -0.013814233 -0.016918911 -0.019536276
# [6] -0.021842220 -0.023926954 -0.025844064 -0.027628466 -0.029304414
# [11] -0.030889564

upper.cone <- S*exp(upper.bounds)
# [1] 23.62000 23.85185 23.94856 24.02303 24.08598 
# [6] 24.14159 24.19197 24.23839 24.28168 24.32241 
# [11] 24.36100

lower.cone <- S*exp(lower.bounds)
# [1] 23.62000 23.39040 23.29595 23.22374 23.16303 
# [6] 23.10968 23.06155 23.01738 22.97635 22.93787 
# [11] 22.90154

foo <- ggplot(data=data.frame(date=index(pcef_ad),data=coredata(pcef_ad)),aes(x=index(pcef_ad),y=coredata(pcef_ad))) + 
  geom_line() + 
  geom_polygon(data=data.frame(date=c(today()+0:n,today()+n:0),data=c(lower.cone,rev(upper.cone))),mapping=aes(x=date,y=data),alpha=1, fill = "pink")

foo
```

### EDIT 2

I tried the code with a "forecast" period of 1000 days, just to see. And there is something in the produced values that I don't understand: Why is the upper range wider than the lower range ? It looks like the prediction takes account of the upside trend but how ? Honestly I expected the range to be symetrical with regard to the last observed price. I expected the blue line and the green line to have the same length, but the green one is bigger ... why ?

## Answer by vanguard2k (score 4, accepted)

https://quant.stackexchange.com/a/37061

First, you are right about your issue with $n$ (the nb_of_days). If you calculate the daily volatility of returns, you only need to multiply it by $n^{1/2}$.

What you are trying to get is the volatility cone. Let $S$ be the last stock price, $p$ the probability of the stock price being below and $Q$ the quantile function of the standard normal.

The distribution of the log of returns is given by:

$$ log(S_1/S_0) \approx N(\mu_s/365,\sigma_S/(365)^{1/2})$$

So, for the whole horizon, the upper border of the cone is given by

$$ S_0 \cdot \text{exp}\left(Q(p)\cdot \sigma_S\cdot (n/365)^{1/2} + \mu_s\cdot n /365\right).$$

The lower bound is the same term with $1-p$ respectively. For shorter time horizons, we just need to change the value for $n$.

It takes some time to produce a good plot but you could try something like this with the ggplot2 package:

```
sigma <- sd(pcef_daily_return)*sqrt(365)
mean <- 0
probability <- 0.95

S <- coredata(tail(pcef_ad,1))[1,1]

upper.bounds <- qnorm(probability,mean=mean*(0:n)/365,sd=sigma*sqrt((0:n)/365))
lower.bounds <- qnorm(1-probability,mean=mean*(0:n)/365,sd=sigma*sqrt((0:n)/365))
upper.cone <- S*exp(upper.bounds)
lower.cone <- S*exp(lower.bounds)

foo <- ggplot(data=data.frame(date=index(pcef),data=coredata(pcef)),aes(x=index(pcef),y=coredata(pcef))) + 
  geom_line() + 
  geom_polygon(data=data.frame(date=c(today()+0:n,today()+n:0),data=c(lower.cone,rev(upper.cone))),mapping=aes(x=date,y=data),alpha=0.25)

foo
```