The jump model we are going to look at is a perfect example of what monte carlo is helpful for. In this model we are going to add in a random jump that happens only once in a while. The jump, when it happens, will draw from a normal distribution but it won’t be scaled by the time increment because it is supposed to be a random large movement that happens with small probability (which we sometimes see in financial markets). The way we do this is to define the jump as:

With the jump dummy, we will define it as 1 if a random uniform number is less than threshold * the time increment. This means that if the threshold is .5, and the time increment is 1/10, then if the random uniform number is less than .05, the variable is 1. Otherwise it is 0. In this context, the threshold standard for the average number of jumps in a year. For this we will say the threshold is .5, which we then will multiply by our time increment of 1/252. This can model it so that on average there is half a jump a year for each run. We can also use a very large sigma of 20% so the jumps can be rather drastic.

Step 1 is to find random uniform numbers.

#For a jump model, the first step is to make a uniform distribution
np.random.seed(1)
dummy = pd.DataFrame(np.random.uniform(0,1,252*10).reshape(252,10))
print(dummy)

Convert to the dummy variable.

#Then we can convert it into dummy variables of 0 or 1 where 1 means there is a jump and 0 means no jump
threshold = .5
time_step = 1/252
dummy = (dummy < (threshold * time_step)).astype(int)
print(dummy)

We can check the frequency of jumps for each of the ten runs and see it happens rarely.

#Note how rare jumps are
print(dummy.sum())

Convert this to a random jump with sigma of 20%.

#Multiply in random jumps
#In the third column, we have a jump of >11% which is represented as a factor of .88 multiplied in
np.random.seed(5)
jump = np.exp(pd.DataFrame(np.random.normal(0,1, 252*10).reshape(252,10)) *.2 - .2**2 /2) * dummy
print(jump)

For adding this to the regular return stream, let's do a few things. First, get rid of anywhere that a jump did not happen by setting it to null. Next, subtract 1 and fill the null values with 0. This is now something we can add on.

#Let's make this into something we can add on
np.random.seed(5)
jump = np.exp(pd.DataFrame(np.random.normal(0,1, 252*10).reshape(252,10)) *.2- .2**2 /2) * dummy
#Set anything without a jump to null then subtract 1 from all jumps
jump = jump.where(jump != 0) - 1
#Fill the null values with 0 because there is no change to the regular returns because there is no jump
jump = jump.fillna(0)
print(jump)

Because our paths variable is not in returns space we need to convert it back to returns space.

#Convert the first 10 paths to returns
pct_returns = paths.iloc[:,:10].pct_change().dropna()
#Reset the index so the jump and regular returns are aligned
pct_returns.index = list(range(len(pct_returns.index)))
print(pct_returns)

Add the two returns.

#Add the two returns
pct_returns = pct_returns + jump

Now plot the paths and see how the jumps add a lot of randomness in. This is important because in real stock markets it tends to happen that crazy events can happen that swing stock prices wildly.

prices = (pct_returns + 1).cumprod()
prices.index = prices.index + 1
prices.loc[0] = 1
prices = prices.sort_index()
prices = prices * 100
prices.plot(kind='line')
plt.xlabel("Time Step")
plt.ylabel("Stock Price")
plt.title("Simulated Paths with Jumps")
plt.show()

Formalize it into a function so we can compare the two!

rf = .02
time_step = 1/252
price0 = 100
sigma = .08
runs = 100000
periods = 252
jump_sigma = .2
jump_prob = .5

def simulate_returns(price0, rf, time_step, periods, sigma):
    returns = pd.DataFrame(np.random.normal(0,1, periods * runs).reshape(periods, runs))
    returns = returns * sigma * timestep ** .5 + (rf - .5 * sigma ** 2) * timestep
    returns = np.exp(returns).cumprod()
    returns.index = returns.index + 1
    returns.loc[0] = 1
    returns = returns.sort_index()
    paths = returns * price0
    return paths
np.random.seed(0)
paths1 = simulate_returns(price0, rf, time_step, periods, sigma)

def simulate_returns_with_jump(price0, rf, time_step, periods, sigma, jump_sigma, jump_prob):
    returns = pd.DataFrame(np.random.normal(0,1, periods * runs).reshape(periods, runs))
    returns = returns * sigma * timestep ** .5 + (rf - .5 * sigma ** 2) * timestep
    
    threshold = jump_prob * time_step
    dummy = pd.DataFrame(np.random.uniform(0,1, periods * runs).reshape(periods, runs))
    dummy = (dummy < threshold).astype(int)
    jump = np.random.normal(0,1, periods * runs).reshape(periods, runs) *jump_sigma- jump_sigma**2 /2
    jump = jump * dummy
    returns = returns + jump
    
    returns = np.exp(returns).cumprod()
    returns.index = returns.index + 1
    returns.loc[0] = 1
    returns = returns.sort_index()
    paths = returns * price0
    return paths
np.random.seed(0)
paths2 = simulate_returns_with_jump(price0, rf, time_step, periods, sigma, jump_sigma, jump_prob)

The mean path stays the same with and without jumps.

pd.concat([paths1.mean(axis=1), paths2.mean(axis=1)], axis=1).plot(kind='line')
plt.legend(['No Jump', 'Jump'])
plt.xlabel("Time Step")
plt.ylabel("Average Stock Price")
plt.title("Average Simulated Paths")
plt.show()

That being said, the volatility is very different. If we find the volatility of log percent change of each run, and then take the mean of that we see this difference. The first set of runs have a volatility of 8% like we specified. The second, however, has a much larger volatility due to these random swings that happen.

vol1 = np.log(paths1.pct_change() + 1).std() * 252 ** .5
print(vol1.mean())
vol2 = np.log(paths2.pct_change() + 1).std() * 252 ** .5
print(vol2.mean())

To drive home this point, look how much longer of a tail the distribution with jumps has. The volatility during some runs can be quite high if the dummy variable is activated often.

import seaborn as sns
sns.kdeplot(vol1)
sns.kdeplot(vol2)

plt.xlabel("Volatility")
plt.ylabel("KDE")
plt.title("Kernel Density Estimation of Volatility")
plt.show()

This will also impact the valuation of call options. Look at what the fair price of a call option might be given a distribution with no jumps and with jumps. Given jumps we notice that a call option should be priced much higher.

print(paths1.iloc[-1].apply(lambda x: max(x-100, 0)).mean() / np.exp(.02))
print(paths2.iloc[-1].apply(lambda x: max(x-100, 0)).mean() / np.exp(.02))

The monte carlo simulations method is quite useful when we want to add unique characteristics to a distribution or we need to apply valuation for a derivative which has a payoff based on what happens not only at the end of the time but also during the time.