# 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


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
# 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

# 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

# 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

# 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"

# 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

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 ?

• Maybe one piece of info from my side. To me, this looks like a plot.xts plot with the areas added - but I dont know why. Try ggplot2! – vanguard2k Nov 24 '17 at 11:11
• Hey @vanguard2k ! Thank you for your interest in my issue. Yes, this plot is just a basic plot made with the plot() function. I added the areas using random values in photoshop to show what I'm looking for because I'm am stuck at calculating those values first. Any idea ? – rmrndr Nov 24 '17 at 11:45
• Hi, its from the xts package and not the standard plot function. The package has an own implementation of the plot function that it takes for objects of type xts. – vanguard2k Nov 24 '17 at 12:57
• Oh ok, I didn't know that – rmrndr Nov 24 '17 at 13:39

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

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

• Thank yo so much I did not know how to address the problem with the qnorm() function. Your code didn't work as is in my console so I made a few adjustment (see my edit in the original post). But thank you so much, you solved all my problems ! – rmrndr Nov 24 '17 at 17:03
• I still have one question though : why do you set a mean at zero ? And why wasn't it needed to convert prices in log somewhere like I did in my exemple ? – rmrndr Nov 24 '17 at 17:05
• I made a second edit at the end of my original post. Do you know why it does that ? – rmrndr Nov 24 '17 at 17:23
• @user_j Well, I set the mean to 0 because I thought thats what you want to do. See, what I am simulating are the log-returns which you assumed to be normal. The mean of their distribution cant be the log of the price since its a return! (guessing (or estimating, as you want) this is the hard part :) ) The cone is not symmetric because you wanted the continuous returns to be normal. After taking the exponential and compounding the stock price, it is not symmetric anymore. Its an exponential of a normal distribution. – vanguard2k Nov 25 '17 at 19:55
• Ok now I was thinking about something else... Why didn't you divide probability by 2 when defining upper.bounds and lower.bounds. The cone defines a range where the price should be 95% Doesn't that mean that it's 47.5% up and 47.5% below? Honestly I don't know and it confuses me... – rmrndr Nov 25 '17 at 21:50