MSF 503 Chapter 8 Simulation 2 PDF

CHAPTER 8 SIMULATION 2: Generating Price Paths A process is a set of steps that converts inputs into outputs. A time series is a sequence of data where each sample or measurement occurs at a consistent time interval. To be sure, financial time series data is not generate...

CHAPTER 8 SIMULATION 2: Generating Price Paths A process is a set of steps that converts inputs into outputs. A time series is a sequence of data where each sample or measurement occurs at a consistent time interval. To be sure, financial time series data is not generated by any process, but we do attempt to model the evolution of financial data using mathematical processes to aid in valuation and prediction. In this chapter, we use the vectorised functionality of the numpy library to generate random numbers concurrently. 8.1 Stochastic Processes A stochastic process is a sequence of random variables. A martingale is a stochastic process where the expected value of xt+1 is equal to xt. For a discrete-time martingale: E ( xt 1 | x1 ,..., xt )  xt (1) A fair coin is one with a probability of heads equal to 0.5. We often assign a numerical value to each flip, such as ( heads, tails ) → ( 1, -1 ). (Notice that, since each sample is a Bernoulli trial, if both the probability of heads and the bet size remain constant, the resulting distribution will be binomial.) The martingale property in (1) ensures that knowledge of past coin flips won’t help predict future ones. This example generates discrete random walks as martingales from a fair coin. import numpy as np import matplotlib.pyplot as plt # Set simulation parameters n_steps = 1000 n_paths = 5 steps = np.zeros( n_paths, n_steps ) # Generate random steps, up or down, with equal probability © 2024 Ben Van Vliet 153  steps = np.random.choice( [ -1, 1 ], size = ( n_paths, n_steps ) ) # Calculate the cumulative sum to get the random walk paths = np.cumsum( steps, axis = 1 ) # Plot the results plt.figure( figsize = ( 12, 6 ) ) for i in range( n_paths ): plt.plot( paths[ i ] ) plt.show() A stochastic process is a Markov process if the conditional probability distribution of future states depends solely upon the present state and is, therefore, independent of all prior states: P( X x1  x | X1  x1 , X 2  x2 ,..., X n  xn )  P( X n1  x | X n  xn ) (2) This is the Markov property, which means the process is “memoryless.” That is, the conditional probability distribution of the next state is dependent only on the current state, and independent of all prior states. (A Markov process need not have an expected future value equal to the current value. Thus, not all martingales satisfy the Markov property and vice versa.) Nearly all financial models satisfy the Markov property. A Markov chain is a stochastic process that transitions from one state to another according to defined transition probabilities. The random transitions satisfy the Markov property because they are functions solely of the transition probabilities of the current state, and not on those of prior states. import numpy as np import matplotlib.pyplot as plt # Define three states and the matrix of transition probabilities states = [ -1, 0, 1 ] © 2024 Ben Van Vliet 154 current_state = 0 transition_matrix = np.array( [ [ 0.7, 0.2, 0.1 ], # Transition probabilities in state 0 [ 0.3, 0.5, 0.2 ], # Transition probabilities in state 1 [ 0.2, 0.3, 0.5 ] # Transition probabilities in state 2 ]) # Initialize the number of steps num_steps = 100 # Start at the current state simulated_states = [ current_state ] for _ in range( num_steps ): # Pick the next state according to the transition probabilities of the current state next_state = np.random.choice( states, p = transition_matrix[ current_state ] ) simulated_states.append( next_state ) current_state = next_state # Calculate the number of times in each state state_counts = [ simulated_states.count( state ) for state in states ] print( "State counts:\n", state_counts ) # Calculate the net probability of being in each state from the simulated data state_probabilities = [ count / num_steps for count in state_counts ] print( "\nSimulated stationary state probabilities:\n", state_probabilities ) # Plot state transitions plt.figure( figsize = ( 12, 6 ) ) plt.plot( simulated_states, marker='o' ) plt.show() © 2024 Ben Van Vliet 155  A diffusion process is a continuous-time Markov process that can generate a continuous random sequence, called a sample path of time series data (such as stock prices). A Levy process is any diffutions process that starts at 0 and has stationary and independent increments. By stationary, we mean that the probability distribution of any segment of the random sequence xt+h − xt depends only upon the length of the time interval h. So, increments with equally long intervals must be identically distributed. The most well-known Levy process is a Wiener process, where xt+h − xt is normally distributed with µ = 0 and σ2 = h, or xt+h − xt ~ N( 0, h ). This variance is important because is shows that the standard deviation is the square root of t, and so it is that volatility scales with the square root of time. We can use a Wiener process to define Brownian motion, where xt equals xt – 1 plus the square root of the time step times a standard normal variate, zs. Wt 1  Wt  t  1  t  z s ,t Clearly here, t + 1 – t = 1, so that as long as the standard deviation is scaled to the correct time interval, the calculation is simplified. That is, one year or one month or one day standard deviation. import numpy as np import matplotlib.pyplot as plt T = 1.0 # Total time is one year N = 250 # Number of days in a year num_paths = 100 # Number of sample paths # Each time step is one day dt = T / N # Wiener process is the standard deviation ( sqrt( time interval ) ) times a random # zs ( i.e. N( mu = 0, var = 1 ) ) dW = np.sqrt( dt ) * np.random.normal( 0, 1, ( num_paths, N ) ) W = np.cumsum( dW, axis = 1 ) W = np.insert( W, 0, 0, axis = 1 ) © 2024 Ben Van Vliet 156  # Plot the results plt.figure( figsize = ( 12, 6 ) ) for i in range( num_paths ): plt.plot( W[ i ] ) plt.show() # Calculate statistics of final values final_values = W[ :, -1 ] print( "\nMean of final values: ", np.mean( final_values ) ) print( "Variance of final values: ", np.var( final_values ) ) print( "Theoretical variance: ", T ) 8.2 Geometric Brownian Motion Geometric Brownian motion (GBM) is a continuous-time stochastic process where the logarithm of the random variable St follows Brownian motion with a drift, according to the following stochastic differential equation: St  S0  dSt   St dt   St dWt (3)  St (  dt   dWt ) where Wt is a Wiener process and μ is the drift rate and σ the volatility. If S is a stock price, this says that the change in the stock price is equal to the price of the stock times a mean drift rate times the change in time plus the standard deviation times some normally distributed random number. The reason we use GBM to model stock price paths is because it encapsulates two widely observed phenomena in financial markets: 1. The long term trends of markets, represented by the mean term. 2. The white noise or random price movements or volatility around the trend, represented by the standard deviation term, which scales with the square root of time. © 2024 Ben Van Vliet 157  Consider compound interest, where number of times the interest compounds per year n → ∞, so that: nt  r e  lim 1   rt (4) n   n Continuous compound interest is: St  Soert (5) The continuous compound interest rate formula is a deterministic process. There is no uncertainty about the future values. But what if we wanted to add some GBM to this process? 8.2.1 Stochastic Differential Equation to Solution Using Ito’s Lemma Given the assumption that stock price follows GBM as in (3). Dividing both sides by the stock price yields: dSt  dt  dWt (6) S0 This shows that the rate of return is equal to the mean drift rate over the change in time plus the stochastic component, standard deviation times the random Wiener value. We can integrate this: t t t dS 0 S0t  0 dt  0 dWt (7) So that: dSt t   (t  0)   (Wt  W0 ) (8) S0 0 And, St  S0  t  Wt (9) Now, let’s start again: f (St )  ln(St ) (10) So, therefore: © 2024 Ben Van Vliet 158  1 1 f ( S t )  and f ( S t )   2 (11) St St Because the graph of the natural log function is non-linear. Now, S  d ln(St )  ln(St )  ln(S0 )  ln t  (12)  S0  We can estimate the change the in the log price of the stock using Taylor’s theorem. We do this precisely because of the non-linearity of the log function. 1 d ln(St )  f ( St )dSt   f ( St )dSt2 (13) 2 Now, Ito’s Lemma states that: lim dSt2  St2 2 dt (14) dt 0 So that by substituting in, we get: 1 d ln(St )  f ( St )dSt   f ( St )  St2 2 dt (15) 2 And furthermore, from equations (11) we get: 1 1 1 d ln(S t )   dSt   2  S t2 2 dt (16) St 2 St Then, 1 1 d ln(S t )   dSt    2 dt (17) St 2 But, from (3) we get: 1 1 d ln(S t )   S t ( dt  dWt )   2 dt (18) St 2 Thus, © 2024 Ben Van Vliet 159  1 d ln(St )  dt  dWt   2 dt (19) 2 And, 1 d ln(St )  (    2 )dt  dWt (20) 2 Via integration, t t t 1 2 0 d ln(St )  0 (  2 )dt  0 dWt (21) So that, t 1 ln(St )|  (    2 )(t  0)   (Wt  W0 ) (22) 0 2 Or, 1 ln(St )  ln(S 0 )  (    2 )t  Wt (23) 2 Then, raising e to the power of each side gives: 1 St (    2 ) t Wt e 2 (24) S0 Solving for St gives: 1 (    2 ) t  Wt St  S 0 e 2 (25) And St is log-normally distributed: 1 St ~ LN ((    2 )t ,  2t ) (26) 2 So that, 1 ln(St ) ~ N ((    2 )t ,  2t ) (27) 2 Thus, ln(St )  ln(S0 ) = μ (28) Now, because the expected value of the Wiener process is 0, the expected value of S is: 1 (    2 ) t  E (Wt ) E ( St )  S 0 e 2 (29) 1 1 (    2 ) t   2t  S0 e 2 2 © 2024 Ben Van Vliet 160   S0et This then shows that using the continuous compounding equation: St  S0ert (30) 8.2.3 Simulating Stock Prices Here is an actual chart using real daily closing price data of Boeing stock (symbol BA) from 10-02-2008 to 10-02-2009. 60 50 40 30 20 10 0 1 13 25 37 49 61 73 85 97 109 121 133 145 157 169 181 193 205 217 229 241 Price path of BA stock The histogram for the daily log returns for BA over the year also appears to be approximately normally distributed: Histogram 45 120.00% 40 100.00% 35 Frequency 30 80.00% 25 60.00% 20 15 40.00% 10 20.00% 5 0 0.00% 1 08 06 04 02 02 04 06 08 1 0. 0. -0 0. 0. 0. 0..... -0 -0 -0 -0 Bin Distribution of BA Log Returns © 2024 Ben Van Vliet 161 From this data we can retrieve the empirical mean and standard deviation. We can generate closing prices using GBM. Let’s start with a single price path. import math import numpy as np import matplotlib.pyplot as plt drift = 0.00 sigma = 0.01 dtime = 1.0 prices = np.zeros( 1000 ) prices[ 0 ] = 100.00 for row in range( 1, 1000 ): prices[ row ] = prices[ row - 1 ] * math.exp( drift * dtime - 0.5 * math.pow( sigma, 2 ) * dtime + sigma * np.random.normal( 0, 1, 1 ) ) plt.plot( prices ) plt.show() Here is the same model using the vectorized functionality of numpy arrays. import numpy as np import matplotlib.pyplot as plt N = 250 # number of time steps dt = 1.0 # time step mu = 0.001 # expected return, or drift, in the correct units sigma = 0.01 # volatility in the correct units S0 = 100.0 # initial price # Generate random increments for the Brownian motion with dt = 1 dW = np.random.normal( 0, 1, N ) © 2024 Ben Van Vliet 162  # Generate the prices using the cumulative sum of the zs S = S0 * np.exp( ( mu - 0.5 * sigma ** 2 ) * dt + sigma * np.cumsum( dW ) ) # Plot the simulated stock price plt.plot( S ) plt.show() Now let’s make 25 price paths. import math import numpy as np import matplotlib.pyplot as plt drift = 0.00 sigma = 0.01 dtime = 1.0 prices = np.zeros( ( 1000, 25 ) ) for i in range( 0, 25 ): prices[ 0 ][ i ] = 100.00 for col in range( 0, 25 ): for row in range( 1, 1000 ): prices[ row ][ col ] = prices[ row - 1 ][ col ] * math.exp( drift * dtime - 0.5 * math.pow( sigma, 2 ) * dtime + sigma * np.random.normal( 0, 1, 1 ) ) plt.plot( prices ) plt.show() © 2024 Ben Van Vliet 163 Here is the same simulation in vectorized code. import numpy as np import matplotlib.pyplot as plt mu = 0.00 # daily drift sigma = 0.01 # daily volatility of returns dt = 1 # Generate daily log returns logret = np.random.normal( mu - 0.5 * sigma ** 2, sigma, size = ( 25, 1000 ) ) # Calculate the price paths price_paths = 100.00 * np.exp( logret.cumsum( axis = 1 ) ) plt.plot( price_paths.T ) plt.show() 8.3 Simulating Short Rates The short rate rt is the instantaneous spot interest rate, expressed as a continuously compounded and annualized number. How does the spot rate evolve over © 2024 Ben Van Vliet 164 time? There are various models. The Vasicek one factor model describes this evolution as being driven by a single factor—Brownian motion—which follows a stochastic, mean- reverting Ornstein-Uhlenbeck process rather than GBM. drt  a(b  rt )dt   dWt Where r0 is the initial interest rate, a is the speed of reversion to the long-run mean b, σ is the instantaneous volatility. The model is often used to value interest rate derivatives. import numpy as np import matplotlib.pyplot as plt r0 = 0.01875 a = 0.2 b = 0.01 sigma = 0.012 T = 10.0 N = 200 dt = T / N rates = [ r0 ] for i in range( N ): dr = a * ( b - rates[ -1 ] ) * dt + sigma * np.random.normal() rates.append( rates[ -1 ] + dr ) plt.plot( range( N + 1 ), rates ) plt.show() © 2024 Ben Van Vliet 165