Basic statistical concepts for Finance
R-bloggers 2020-05-24
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Returns VS Log-Returns
Load the quantmod package to retrieve financial data (see tutorial on Data Acquisition in R).
Get Apple stock quotes (daily).
Returns
\[R_t = \frac{P_t-P_{t-1}}{P_{t-1}} = \frac{P_t}{P_{t-1}}-1\]
Note the use of the lag function in the following code snippet. The function shifts the time base of a xts object back by a given number of observations. Warning: the function behaves differently for zoo, numeric or other object classes.
# shift the series back by one period and compute returnsr <- AAPL$AAPL.Close/lag(AAPL$AAPL.Close, k = 1) - 1Log-Returns
\[r_t = ln(P_t)-ln(P_{t-1})\]
# shift the series back by one period and compute log-returnsr.log <- log(AAPL$AAPL.Close) - log(lag(AAPL$AAPL.Close, k = 1))Returns cannot be computed for the first observation, as there is no previous observation to use. R returns NA in this case. Clean the data.
Log-returns are approximately equal to standard returns when they are small. Note: approximately does not mean equal. \[r_t=ln(P_t)-ln(P_{t-1})=ln\Bigl(\frac{P_t}{P_{t-1}}\Bigl)=ln\Bigl(\frac{P_t-P_{t-1}+P_{t-1}}{P_{t-1}}\Bigl)=ln(1+R_t) \approx R_t\] where the last passage can be easily obtained by Taylor expansion. Compare the returns.
# build a data.frame containing ret and log-retx <- data.frame(r, r.log) colnames(x) <- c('ret', 'log.ret') head(x) ## ret log.ret## 2007-01-04 0.022195679 0.021952941## 2007-01-05 -0.007121188 -0.007146665## 2007-01-08 0.004938272 0.004926118## 2007-01-09 0.083070106 0.079799700## 2007-01-10 0.047855589 0.046745780## 2007-01-11 -0.012371021 -0.012448179Plot the returns. For high quality charts, refer to the ggplot2 package.
# base R plot functionplot(x$log.ret ~ x$ret, main = 'Returns VS Log-Returns', ylab = 'Log-Returns', xlab = 'Returns', xlim = c(-0.21, 0.21), ylim = c(-0.21, 0.21))# add the line y = x abline(b = 1, a = 0, col = 'green')
Multi-Period Returns
\[R_{t,s}=\frac{P_t}{P_s}-1=\frac{P_t}{P_{t-1}}\frac{P_{t-1}}{P_s}-1=\frac{P_t}{P_{t-1}}\frac{P_{t-1}}{P_{t-2}}…\frac{P_{s+1}}{P_s}-1=\\ =(R_{t}+1)(R_{t-1}+1)…(R_{s+1}+1) – 1=\\ =\Bigl(\prod_{i=s+1}^{t}(R_i+1)\Bigl)-1\]
## [1] 0.498575# compound calculationr.i <- window(r, start = '2016-01-01', end = '2018-12-31') # extract returnsprod(r.i+1) - 1## [1] 0.498575Multi-Period Log-Returns
\[r_{t,s}=ln(P_t)-ln(P_s)=ln(P_t)-ln(P_s)=ln(P_t)-ln(P_{t-1})+ln(P_{t-1})-ln(P_s)=\\ =ln(P_t)-ln(P_{t-1})+ln(P_{t-1})-ln(P_{t-2}) … +\;ln(P_{s+1})-ln(P_s)=\\ =r_t+r_{t-1}+…\;+r_{s+1}=\sum_{i=s+1}^tr_i\]
# direct calculation log(AAPL$AAPL.Close['2018-12-31'][[1]]) - log(AAPL$AAPL.Close['2015-12-31'][[1]])## [1] 0.4045146# compound calculationr.i <- window(r.log, start = '2016-01-01', end = '2018-12-31') # extract returnssum(r.i)## [1] 0.4045146Cumulative Returns
Compute the standard returns with respect to the first observation at each point in time (\(R_{t,0}\) for every \(t\)). Use the following approaches:
- (RIGHT) cumulate standard daily returns
- (RIGHT) cumulate log-returns and convert to standard return using the following property: \(r=ln(1+R)\rightarrow R=exp(r)-1\)
- (WRONG) cumulate log-returns as if they were standard returns. This holds approximately, and the error becomes larger and larger when cumulating more and more returns
# cumulate standard returnsr.cum <- cumprod(1+r) - 1# cumulate log-returns and convert to standard returnsr.log.cum.right <- exp(cumsum(r.log)) - 1# cumulate log-returns as if they were standard returnsr.log.cum.wrong <- cumprod(1+r.log) -1 # create a unique xts object x <- merge(r.cum, r.log.cum.wrong, r.log.cum.right)colnames(x) <- c('R', 'r.wrong', 'r.right')# plot and compareplot(x, legend.loc = 'topleft', main = 'Cumulative return computation')
The green line coincides with the black one. This can be checked printing the values.
## R r.wrong r.right## 2007-01-04 0.02219568 0.02195294 0.02219568## 2007-01-05 0.01491643 0.01464939 0.01491643## 2007-01-08 0.01992836 0.01964767 0.01992836## 2007-01-09 0.10465392 0.10101525 0.10465392## 2007-01-10 0.15751779 0.15248306 0.15751779## 2007-01-11 0.14319811 0.13813675 0.14319811Adjusted Price
The closing price is the ‘raw’ price which is just the cash value of the last transacted price before the market closes. Adjusted closing price amends a stock’s closing price to accurately reflect that stock’s value after accounting for any corporate actions (e.g. dividends, splits). It is considered to be the true price of that stock and is often used when examining historical returns or performing a detailed analysis of historical returns. Compute returns using the adjusted price.
# compute returns using adjusted pricesr.adj <- AAPL$AAPL.Adjusted/lag(AAPL$AAPL.Adjusted, k = 1) - 1# drop NAr.adj <- na.omit(r.adj)# cumulate returnsr.adj.cum <- cumprod(1+r.adj) - 1# add adjusted returns to the ones computed in the previous sectionx <- merge(x, r.adj.cum)colnames(x)[4] <- 'R.adj'# plot and compareplot(x, legend.loc = 'topleft', main = 'Cumulative return computation')
Central Limit Theorem
The central limit theorem establishes that when independent and identically distributed (i.i.d.) random variables are added, their properly normalized sum tends toward a normal distribution even if the original variables themselves are not normally distributed.
\[\frac{S_n-\mu_n}{\sigma_n}\rightarrow N(0,1)\] Define the following function to test the Central Limit Theorem
# n: number of iid random variables to sum up# plot: plotting results? Default FALSECLT <- function(n, plot = FALSE){ # number of trials N <- 100000 # build a matrix with # - 'N' rows (numer of trials) # - 'n' columns (number of iid random variables to sum up) # - 'n*N' uniform random variables: 'runif' function u <- matrix(data = runif(n = n * N), ncol = n, nrow = N) # for each row, sum up all columns -> generate Sn z <- rowSums(u) # normalize according to the theorem z <- (z - mean(z))/sd(z) # if plot if(plot){ # generate histogram of Sn hist(z, freq = FALSE, breaks = 100) # add normal distribution x <- seq(-3,3,0.01) lines(x = x, y = dnorm(x = x), col = 'blue') } # else return Sn values else{ return(z) } }Run the function with several values for n. The distribution of Sn converges to a Normal distribution as n increases.
Normality of financial returns
Plot AAPL log-returns. Are they Normal?
# histogram of AAPL log-returnshist(r.log, breaks = 50, main = "Density of Log-Returns", xlab = 'Log-Returns', freq = FALSE)# mean r.log.mean <- mean(r.log)# standard deviationr.log.sd <- sd(r.log) # density of N(mu, sigma)x <- seq(min(r.log), max(r.log), by = 0.001) y <- dnorm(x = x, mean = r.log.mean, sd = r.log.sd) # add fitted normal distributionlines(x = x, y = y, col = 'blue')
AAPL daily log-returns, and stock returns in general, are not Normal (Fama 1965, Journal of Business, Nobel price 2013). On the other hand, daily log-returns can be thought as the cumulative sum of intraday log-returns and, according to the central limit theorem, they should be normally distributed. This means that intraday stock returns are not iid and the central limit theorem does not hold in this case.
Law of Large Numbers
The law of large numbers is a theorem that describes the result of performing the same experiment a large number of times. According to the law, the average of the results obtained from a large number of trials should be close to the expected value, and will tend to become closer as more trials are performed. In other words, the mean is a good estimator of the expected value.
\[\bar{X_n} \rightarrow E[X]\] Estimate the expected value of a uniform random variable that takes values in [0,1].
## [1] 0.4919098# generate 100000 uniform random variables to increase precision mean(runif(n = 100000, min = 0, max = 1))## [1] 0.5003542For simulated data it is possible to increase the precision of the estimation by increasing the number of trials. For real data, where the sample size is finite, the precision cannot be increased arbitrarily and we need to understand how uncertain the estimation is. In other words, we need an estimation of the precision, together with the estimation of the mean. This is usually called standard error or standard deviation of the mean. When realized returns are iid, the standard error of the mean can be easily computed by \(\sigma_n/\sqrt{n}\), where \(\sigma_n\) is the standard deviation of the realized returns and \(n\) is the sample size. When realized returns are not iid, it is under certain conditions (“stationarity of returns”) still possible to compute the standard error, but the required formula is slightly more complicated. Moreover, the estimation of the mean, for large \(n\), is normally distributed according to the central limit theorem.
\[\bar{X_n} \rightarrow N\Bigl(\bar{X_n}, \frac{\sigma_n}{\sqrt{n}}\Bigl)\] Estimate the expected daily return for AAPL and its distribution.
# estimatation of the expected daily return mu <- mean(r)# estimation of the standard deviation of the samplesigma <- sd(r)# sample sizen <- length(r)# standard devation of the meanst.err <- sigma/sqrt(n)# plot the estimated distribution of the expcted return x <- seq(mu-5*st.err, mu+5*st.err, length.out = 1000)plot(x = x, y = dnorm(x = x, mean = mu, sd = st.err), type = 'l', ylab = 'Density', xlab = 'Expected Daily Return')
## [1] 0.001181849# how likely is the expected return to be less than 0.05%?pnorm(q = 0.0005, mean = mu, sd = st.err, lower.tail = TRUE)## [1] 0.02580614Remark. The law of large numbers holds when performing the same experiment a large number of times. This is not the case for financial and economic data. Assuming returns coming from the same experiment every day ignores the different economic conditions, sentiment and environment of the different periods.
Stochastic Process
A stochastic process can be defined as a collection of random variables.
Brownian Motion
Increments are normally distributed:
\[\Delta X_t \sim N(\mu, \sigma)\]
The following code snippet defines a function to simulate trajectories of a Brownian Motion.
# n: number of trajectories to simulate# t: number of periods for a single trajectory# mu: mean of the increments# sigma: standard deviation of the increments# x.0: initial value# plot: plot trajectories? default TRUEBM.sim <- function(n, t = 252, mu = 0.001, sigma = 0.05, X.0 = 1, plot = TRUE){ # generate increments and store in a matrix with: # - 't' rows (numer of time points for each trajectory) # - 'n' columns (number of trajectories to simulate) r <- matrix(data = rnorm(n = n*t, mean = mu, sd = sigma), ncol = n, nrow = t) # for each column (trajectory), sum increments starting from x.0 X.t <- apply(r, MARGIN = 2, function(x) X.0 + cumsum(x)) # plot and return if(plot) plot(as.zoo(X.t), screens = 'single', ylab = 'Trajectory') return(X.t)}Simulate the Brownain Motion.
The Brownian Motion can assume negative values. This is not well-suited for financial applications.
Geometric Brownian Motion
Returns are normally distributed. More precisely, log returns are normally distributed and simple returns are log-normal distributed:
\[\frac{\Delta X_t}{X_t} \sim N(\mu, \sigma)\] The following code snippet defines a function to simulate trajectories of a Geometric Brownian Motion.
# n: number of trajectories to simulate# t: number of periods for a single trajectory# mu: mean of the returns# sigma: standard deviation of the returns# x.0: initial value# plot: plot trajectories? default TRUE# log.scale: plot on a log-scale? default FALSEGBM.sim <- function(n, t = 252, mu = 0.001, sigma = 0.05, X.0 = 1, log.scale = FALSE, plot = TRUE){ # generate returns and store in a matrix with: # - 't' rows (numer of time points for each trajectory) # - 'n' columns (number of trajectories to simulate) r <- matrix(data = rnorm(n = n*t, mean = mu, sd = sigma), ncol = n, nrow = t) # for each column (trajectory), cumulate returns starting from x.0 X.t <- apply(r, MARGIN = 2, function(x) X.0 * cumprod(1+x)) # plot and return if(plot){ if(log.scale) plot(as.zoo(X.t), screens = 'single', log = 'y', ylab = 'Trajectory') else plot(as.zoo(X.t), screens = 'single', ylab = 'Trajectory') } return(X.t)}Simulate the Geometric Brownain Motion.
The Geometric Brownian Motion is positive defined. Plot the trajectories on a log-scale and check.
# set the same RNG seedset.seed(123)# simulate 100 trajectories and plot on a log-scalegbm <- GBM.sim(100, log.scale = TRUE)
Exercise
Compute the mean, standard deviation, skewness, and kurtosis for both returns and log-returns for each of the seven investment strategies that can be found in this file. The file contains stock data from the website of Kenneth R. French. Which is your favourite strategy?
# import datadata <- read.csv('https://storage.guidotti.dev/course/asset-pricing-unine-2019-2020/basic-statistical-concepts-for-finance.csv', sep = ';', stringsAsFactors = FALSE)# convert data to xts data <- xts(data[,-1], order.by = as.yearmon(data[,1]))# returns and log-returnsr <- data/lag(data,1)-1 r.log <- log(data)-log(lag(data,1))# drop NAr <- na.omit(r)r.log <- na.omit(r.log)# skewness function skewness <- function(x){ mean((x-mean(x))^3)/sd(x)^3}# kurtosis functionkurtosis <- function(x){ mean((x-mean(x))^4)/var(x)^2 - 3}## Market SizeLo30 SizeMed40 SizeHi30 ValueLo30 ValueMed40 ValueHi30 ## 0.1171537 0.2089651 0.1490341 0.1239627 0.1138287 0.1238843 0.1592684## Market SizeLo30 SizeMed40 SizeHi30 ValueLo30 ValueMed40 ValueHi30 ## 0.09315030 0.13532657 0.10865334 0.09704191 0.09021366 0.09727667 0.11893181## Market SizeLo30 SizeMed40 SizeHi30 ValueLo30 ValueMed40 ValueHi30 ## 0.2003840 0.4149276 0.2766021 0.2131912 0.2012226 0.2102552 0.2709719## Market SizeLo30 SizeMed40 SizeHi30 ValueLo30 ValueMed40 ValueHi30 ## 0.1946990 0.3334365 0.2536980 0.2071972 0.1933582 0.2078262 0.2490666## Market SizeLo30 SizeMed40 SizeHi30 ValueLo30 ValueMed40 ValueHi30 ## -0.39025418 1.40580597 0.38599846 -0.27243948 -0.31865990 -0.08623006 0.38710970## Market SizeLo30 SizeMed40 SizeHi30 ValueLo30 ValueMed40 ValueHi30 ## -0.9599914 -0.1495494 -0.6381693 -1.0442356 -0.7637205 -1.4922693 -0.8256848## Market SizeLo30 SizeMed40 SizeHi30 ValueLo30 ValueMed40 ValueHi30 ## -0.08238205 4.81097118 1.53081600 0.44629946 -0.49069275 2.11458192 1.94780192## Market SizeLo30 SizeMed40 SizeHi30 ValueLo30 ValueMed40 ValueHi30 ## 1.0559344 0.5799876 0.8924158 1.7487555 0.2067589 5.5715657 1.9045128# cumulative returns for each strategyr.cum <- apply(r, MARGIN = 2, function(x) cumprod(1+x))# convert to xts r.cum <- xts(r.cum, as.yearmon(rownames(r.cum)))# plotplot(r.cum, legend.loc = 'topleft')
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.