Assignment 7

Objective and Overview

Objective: Learn to use the numpy numerical programming toolkit and Monte Carlo simulation. Additional practice with object-oriented programming, class definintions, inheritance, and polymorphism.

Overview: In this assignment, you will implement some classes to simulate stock returns and price movements (task 1) to assist with pricing path-dependent options (task 2). In the first task, you will implement a base class to encapsulate the fundamental data members and common calculations used to simulate stock returns using Monte Carlo simulation. In the second task, you will create subclasses to implement several option pricing algorithms.

Preliminaries

In your work on this assignment, make sure to abide by the collaboration policies of the course.

If you have questions while working on this assignment, please post them on Piazza! This is the best way to get a quick response from your classmates and the course staff.

pip3 install numpy scipy matplotlib pandas

Task 1: Simulating Stock Returns

50 points; individual-only

Monte Carlo Simulation is a mathematical technique in which we use data generated at random to help understand what might occur in some real scenario. It is used in a wide variety of domains to simulate outcomes that follow some known (or assumed) probability distribution. In finance, Monte Carlo methods are used to simulate stock returns, and can be instrumental in pricing some assets for which an analytical pricing formula does not exist.

In this task, you will write a definition for the class MCStockSimulator, which encapsulates the data and methods required to simulate stock returns and values. This class will serve as a base class for option pricing classes (in part 2).

Do this entire task in a file called a7task1.py.

1. Create the __init__ and __repr__ methods for the class MCStockSimulator, such that you can create and print out an object like this:

   Example:

   >>> # initial stock price = $100; 1 year time frame
   >>> # expected rate of return = 10%; standard deviation = 30%; 
   >>> # 250 discrete time periods per year
   >>> sim = MCStockSimulator(100, 1, 0.1, 0.3, 250)
   >>> print(sim)
   StockSimulator (s=$100.00, t=1.00 (years), r=0.10, sigma=0.30, nper_per_year=250
   

   The data required include:
   s (the current stock price in dollars),
   t (the option maturity time in years),
   r (the annualized rate of return on this stock),
   sigma (the annualized standard deviation of returns),
   nper_per_year (the number of discrete time periods per year)
   

2. Write a method generate_simulated_stock_returns(self) on your class MCStockOption, which will generate and return a np.array (numpy array) containing a sequence of simulated stock returns over the time period t.

   For example, here is a MCStockSimulator for a 1-year time horizon, with 2 discrete periods per year:

   >>> sim = MCStockSimulator(100, 1, 0.10, 0.30, 2)
   >>> sim.generate_simulated_stock_returns()
   array([ 0.10370361,  0.04763965])
   

   Here is a MCStockSimulator for a 1-year time horizon, with 4 discrete periods per year:

   >>> sim = MCStockSimulator(100, 1, 0.10, 0.30, 4)
   >>> sim.generate_simulated_stock_returns()
   array([ 0.06111529, -0.03161354,  0.01123988, -0.02747573])
   

   Here is a MCStockSimulator for a 0.5-year time horizon, with 2 discrete periods per year:

   >>> sim = MCStockSimulator(100, 0.5, 0.10, 0.30, 2)
   >>> sim.generate_simulated_stock_returns()
   array([ 0.01746703])
   

   Notice that in 0.5 years, it only generates one simulated return (i.e, 2 periods per for 0.5 years is one periodic return).

   Finally, here is a 5-year simulation with 250 discrete time periods per year (i.e., 250 stock market trading days per year).

   >>> sim = MCStockSimulator(100, 5, 0.10, 0.30, 250)
   >>> returns = sim.generate_simulated_stock_returns()
   >>> len(returns)
   1250
   

   Algorithm to generate simulated stock returns

   One method for creating simulated stock returns is to assume that future returns will follow the same probability distribution as historical returns. Let’s assume that stock returns have historically been approximately normally distributed with a mean mu and standard deviation sigma. We can simiulate the annual rate of return as:

   

   where Z is a randomly-drawn number from the standard normal distribution. You can draw such a number using the function np.random.normal().

   To find a simulated rate of return for an arbitrary discrete time period dt, we can use this formula:

   

   Where mu is the mean historical rate of return on the stock, and sigma is the historical standard deviation on the stock.

   Notes/Hints:

   • Each run of generate_simulated_stock_returns() will return a new, independent sequence of simulated stock returns. You should not expect to get the same result twice!

   • The returns in this sequence represent a simluated run through t years, with nper_per_year discrete steps per year. It might be helpful to create a variable dt = 1/nper_per_year to represent the length of each discrete time period.

3. Write a method generate_simulated_stock_values(self) on your class MCStockOption, which will generate and return a np.array (numpy array) containing a sequence of stock values, corresponding to a random sequence of stock return.

   The sequence of stock values is created as follows:

   

   We refer to this sequence of prices as the price path of the stock, i.e., the price of the stock at each discrete time period from now until the end of the simulation. There are t * nper_per_year discrete time periods.

   For example, here is a 1-year simulation with the stock value in each of 4 discrete time periods:

   >>> sim = MCStockSimulator(100, 1, 0.10, 0.30, 4)
   >>> sim.generate_simulated_stock_values()
   array([ 100.,  106.21081853,   99.06516915,  106.82487454, 111.7185901 ])
   

   For example, here is a 2-year simulation with the stock value in each of 24 discrete time periods:

   >>> sim = MCStockSimulator(100, 2, 0.10, 0.30, 24)
   >>> sim.generate_simulated_stock_values()
   array([ 100. ,   98.43464757,   89.06787115,   89.45335537,
    98.10990421,   95.47948379,   98.58499305,  102.60564175,
   110.89048152,  109.14211213,  101.505116  ,  102.75341162,
    94.28467862,   95.64829725,  106.66955089,  116.07599088,
   108.20518856,  103.25591899,  117.17127503,  101.43353659,
    97.12297512,  100.56849565,   98.02435776,  103.12016945,
   105.77297333,  110.350575  ,  112.5208508 ,  121.42911882,
   122.11935187,  125.16098769,  129.47382744,  149.55685873,
   147.32158875,  153.70553074,  160.91099803,  166.92453169,
   161.94170337,  159.82663272,  158.93655807,  144.25364641,
   137.70097581,  140.1791159 ,  139.74213887,  143.07955822,
   134.06674212,  127.83268835,  121.31988325,  131.96440002,
   136.10952943])
   

   Notes/Hints:

   • You should re-use your generate_simulated_stock_returns() to obtain the sequence of returns for the simulation. Next, you will use the returns to create the stock values for the simulation, beginning with the initial stock price s.

4. Write the method plot_simulated_stock_values(self, num_trials = 1), that will generate a plot of of num_trials series of simulated stock returns. num_trials is an optional parameter; if it is not supplied, the default value of 1 will be used.

   For example, consider this sequence of operations:

   >>> sim = MCStockSimulator(100, 2, 0.10, 0.30, 24)
   >>> sim.plot_simulated_stock_values()
   

   which generates this plot:

   

   We can also generate multiple simulations on a single plot, by passing in the parameter num_trials. For example:

   >>> sim = MCStockSimulator(100, 2, 0.10, 0.30, 250)
   >>> sim.plot_simulated_stock_values(5)
   

   which generates this plot:

   

   In the above plot, each line represents one simulated price-path over a period of 5 years, with 250 discrete time periods per year.

   Notes/Hints:

   • Each run of plot_simulated_stock_values() will return a new, independent plot of simulated stock values. You should not expect to get the same result twice!

   • You should re-use your generate_simulated_stock_values() method to obtain a sequence of stock values for each of num_trials trials.

   • Use the matplotlib library to generate the plots.

   • The matplotlib plot method accepts parameters that are of type numpy.array.
     We can pass in a 2-dimensional numpy.array of y-values, where each y-value corresponds to one x-value.

   • The matplotlib library will automatically assign colors to each plot line, so you don’t need to.

Task 2: Pricing Path-Dependent Options

50 points; individual-only

The Black-Scholes method provides the most efficient way to value standard options, such as the European call or put options. However, there exist many exotic options types for which the payoff depends on the price path of the underlying assets, and thus require a technique to simulate the path of a stock price (or other underlying asset) through time.

We will use object-oriented techniques to implement a hierarchy of option classes, extensible to many different kinds of optoins with similar fundamental characteristics such as the underlying stock’s mean rate of return and standard deviation of returns, but with different payoff algorithms (e.g., European, average price, or no-regret options).

Do this task in a file called a7task2.py. You will need to import the class MCStockSimulator from your a7task1.py file, as well as importing the numpy module, i.e.,

    from a7task1 import MCStockSimulator
    import numpy as np

1. Write a class definition for the class MCStockOption, which inherits from MCStockSimulator. This class will encapsulate the idea of a Monte Carlo stock option, and will contain some additional data members( that are not part of class MCStockSimulator and are required to run stock-price simulations and calculate the option’s payoff. However, this class will not be a concrete option type, and thus will not implement the option’s value algorithm, which differs for each concrete option type.

   a. Write an __init__ method that takes the following parameters:
   s, which is the initial stock price
   x, which is the option’s exercise price
   t, which is the time to maturity (in years) for the option
   r, which is the (expected) mean annual rate of return on the underlying stock
   sigma, which is the annual standard deviation of returns on the underlying stock
   nper_per_year, which is the number of discrete time periods per year with which to evaluate the option, and
   num_trials, which is the number of trials to run when calculating the value of this option
   

   The constructor for MCStockOption will need to explicitly call the super class’s constructor to initialize the data attributes of the super class. Only two new data attributes are created in this subclass.

   b. Write a __repr__ method to create a nicely formatted printout of the MCStockOption object, which will be useful for debugging your work.

   >>> option = MCStockOption(90, 100, 1.0, 0.1, 0.3, 250, 10)
   >>> print(option)
   MCStockOption, s=90.00, x=100.00, t=1.00, r=0.10, sigma=0.30, nper_per_year=250, num_trials=10
   

2. Add the following methods to your class MCStockOption:

   a. the method value(self), which will return the value of the option. This method cannot be concretely implemented in this base class, but will be overridden in each subclass (see below). For now, the methods should print out a message and return 0.

   Example:

   >>> option = MCStockOption(90, 100, 1.0, 0.1, 0.3, 250, 10)
   >>> option.value()
   Base class MCStockOption has no concrete implementation of .value(). # print statement
   0 # this is the return value
   

   b. the method stderr(self), which will return the standard error of this option’s value. The standard error is calculated as stdev / sqrt(num_trials), where stdev is the standard deviation of the values obtained from many trials.

   The standard error can only be calculated after running some trials, and as such it cannot be calculated until we implement the value(self) method in the concrete subclasses (explained below). For now, you should use the following method definition in class MCStockOption:

   def stderr(self):
       if 'stdev' in dir(self):
           return self.stdev / math.sqrt(self.num_trials)
   return 0
   

3. Write a class definition for the class MCEuroCallOption which inherits from the base class MCStockOption.

   For example:

   >>> call = MCEuroCallOption(90, 100, 1, 0.1, 0.3, 100, 1000)
   >>> print(call)
   MCEuroCallOption, s=90.00, x=100.00, t=1.00, r=0.10, sigma=0.30, nper_per_year=100, num_trials=1000
   >>> call.value()
   10.164874441384734 ## your value WILL be different but should be relatively close
   

   Calculating the option value:

   The value of the European call option is calculated by:

   

   where S_t is the value of the underlying stock in the last discrete time period.

   Note that each time your execute the .value() method, you will get a different result to to the nature of the random trials:

   >>> call.value()
   11.322180528696247
   >>> call.value()
   10.748354948112402
   >>> call.value()
   10.695278628009447
   

   The standard error of the option’s value describes how good or bad this estimate is. In this case for a call option value calculated with 1000 random trials, the standard error is quite large:

   >>> call.stderr()
   0.5787576231663906
   

   We can obtain a better estimate of the option value (i.e., with lower standard error) by increasing the number of trials:

   >>> call = MCEuroCallOption(90, 100, 1.0, 0.1, 0.3, 100, 1000)
   >>> print(call)
   MCEuroCallOption, s=90.00, x=100.00, t=1.00, r=0.10, sigma=0.30, nper_per_year=100, num_trials=1000
   >>> call.value()
   11.000330304756181
   >>> call.stderr()
   0.6023949256312296
   >>> # change num_trials from 1,000 to 100,000
   >>> call.num_trials = 100000
   >>> call.value() # note: this took about 8 seconds to run on my computer
   10.512108280184108
   >>> call.stderr()
   0.05937301178540991
   

   The value of this option given by Black-Scholes is:

   >>> bs_call = BSEuroCallOption(100, 100, 1.0, 0.3, 0.10)
   >>> bs_call.value()
   10.51985812604488  # the MC option value was pretty close
   

   We can obtain a closer estimate by further increasing the number of trials, i.e.,

   >>> # change num_trials from 1,000 to 1,000,000
   >>> call.num_trials = 1000000
   >>> call.value() # this calculation took about 1 minutes on my computer.
   10.537749946259456
   >>> call.stderr()
   0.01878190454494842
   

   Notice that using Monte Carlo simulation to calculate the value is very inefficient for European call/put options. However, pricing European call options via our Monte Carlo simulation is a good way to check that our simulation is working correctly, i.e., we can check the value obtained for these options against the values from Black-Scholes.

   Here is one more example against which you can test your work:

   >>> bs_call = BSEuroCallOption(40, 40, 0.25, 0.3, 0.08)
   >>> bs_call.value()
   2.7847366578216608
   >>> call = MCEuroCallOption(40, 40, 0.25, 0.08, 0.30, 100, 100000)
   >>> print(call)
   MCEuroCallOption, s=40.00, x=40.00, t=0.25, r=0.08, sigma=0.30, nper_per_year=100, num_trials=100000
   >>> call.value()
   2.7946173763090427
   >>> call.stderr()
   0.012950161471655414
   

4. Write a class definition for the class MCEuroPutOption which inherits from the base class MCStockOption.

   For example:

   >>> put = MCEuroPutOption(100, 100, 1.0, 0.1, 0.3, 100, 1000)
   >>> print(put)
   MCEuroPutOption, s=100.00, x=100.00, t=1.00, r=0.10, sigma=0.30, nper_per_year=100, num_trials=1000
   >>> put.value()
   7.2857299607715076
   

   Calculating the option value:

   The value of the European put option is calculated by:

   where S_t is the value of the underlying stock in the last discrete time period.

   Note that each time your execute the .value() method, you will get a different result to to the nature of the random trials:

   >>> put.value()
   7.1295404189273102
   >>> put.value()
   6.8666625992759549
   >>> put.value()
   7.140392750931456
   >>> put.value()
   7.1207833023244609
   

   The standard error of the option’s value describes how good or bad this estimate is. In this case for a call option value calculated with 1000 random trials, the standard error is quite large:

   >>> put.stderr()
   0.360437099279606
   

   We can obtain a better estimate of the option value (i.e., with lower standard error) by increasing the number of trials:

   >>> put = MCEuroPutOption(100, 100, 1.0, 0.1, 0.3, 100, 100000)
   >>> print(put)
   MCEuroPutOption, s=100.00, x=100.00, t=1.00, r=0.10, sigma=0.30, nper_per_year=100, num_trials=100000
   >>> put.value() ## this took about 8 seconds on my computer
   7.2195765228911499
   >>> put.stderr()
   0.03564371182285262
   

   The value of this option given by Black-Scholes is:

   >>> bs_put = BSEuroPutOption(100, 100, 1.0, 0.3, 0.10)
   >>> bs_put.value()
   7.2178753859826088
   

   Again, we note that Monte Carlo simulation to calculate the value is very inefficient for European call/put options. However, pricing European put options via our Monte Carlo simulation is a good way to check that our simulation is working correctly, i.e., we can check the values obtained for these options against the values from Black-Scholes.

5. Write a class definition for the class MCAsianCallOption which inherits from the base class MCStockOption. The Asian call option’s payoff the amount by which the average stock price exceeds the option’s exercise price during the option’s lifetime.

   (The “Asian” or average price option was developed by two British bankers, who just happened to be in Tokyo when “they developed the first commercially used pricing formula for options linked to the average price of crude oil.” - Wikipedia)

   Calculating the option value:

   The value of the Asian call option (in a single trial) is calculated by:

   where S is price-path of the of the underlying stock form now until time t.

   For example:

   >>> acall = MCAsianCallOption(100, 100, 1.0, 0.10, 0.30, 100, 1000)
   >>> print(acall)
   MCAsianCallOption, s=100.00, x=100.00, t=1.00, r=0.10, sigma=0.30, nper_per_year=100, num_trials=1000
   >>> acall.value()
   8.7286677639144141
   >>> acall.stderr()
   0.37251965838780865
   

   We can obtain a better estimate of the option value (i.e., with lower standard error) by increasing the number of trials:

   >>> # change num_trials to reduce the standard error:
   >>> acall.num_trials = 100000
   >>> acall.value()
   9.082241249588852
   >>> acall.stderr()
   0.039660682189637385
   

   There is no “exact” right answer, but we can increase the number of trials until the standard error is sufficiently small (e.g., maybe a standard error of $0.01 is “good enough”.)

   Here are some additional examples of Asian call options with different parameters:

   >>> acall = MCAsianCallOption(35, 30, 1.0, 0.08, 0.25, 100, 100000)
   >>> acall.value()
   6.10676499229242
   >>> acall.stderr()
   0.014791094213907266
   >>> acall = MCAsianCallOption(35, 40, 1.0, 0.08, 0.40, 100, 100000)
   >>> acall.value()
   1.9336597090613712
   >>> acall.stderr()
   0.013909777807322249
   

6. Write a class definition for the class MCAsianPutOption which inherits from the base class MCStockOption. The Asian put option’s payoff the amount by which the option’s exercise price exceeds average stock price during the option’s lifetime.

   Calculating the option value:

   The value of the Asian put option (in a single trial) is calculated by:

   where S is price-path of the of the underlying stock form now until time t.

   For example:

   >>> aput = MCAsianPutOption(100, 100, 1.0, 0.10, 0.30, 100, 1000)
   >>> print(aput)
   MCAsianPutOption, s=100.00, x=100.00, t=1.00, r=0.10, sigma=0.30, nper_per_year=100, num_trials=1000
   >>> aput.value()
   4.3344210537230685
   >>> aput.stderr()
   0.2126994125989547
   

   We can obtain a better estimate of the option value (i.e., with lower standard error) by increasing the number of trials:

   >>> # change num_trials to reduce the standard error:
   >>> aput.num_trials = 100000
   >>> aput.value()
   4.3828743625359667
   >>> aput.stderr()
   0.022065988231296036
   

   There is no “exact” correct value, but we can increase the number of trials until the standard error is sufficiently small (e.g., maybe a standard error of $0.01 is “good enough”.)

   Here are some additional examples of Asian put options with different parameters:

   >>> aput = MCAsianPutOption(35, 30, 1.0, 0.08, 0.25, 100, 100000)
   >>> aput.value()
   0.17081451092121083
   >>> aput.stderr()
   0.0021427462876315935
   >>> aput = MCAsianPutOption(35, 40, 1.0, 0.08, 0.40, 100, 100000)
   >>> aput.value()
   5.202096070272916
   >>> aput.stderr()
   0.01579490624224986
   

7. Write a class definition for the class MCLookbackCallOption which inherits from the base class MCStockOption. The look-back call option’s payoff the amount by which the option’s exercise price exceeds maximum stock price during the option’s lifetime.

   Calculating the option value:

   The value of the look-back call option (in a single trial) is calculated by:

   where S is price-path of the of the underlying stock form now until time t.

   For example:

   >>> lcall = MCLookbackCallOption(100, 100, 1.0, 0.10, 0.30, 100, 1000)
   >>> print(lcall)
   MCLookbackCallOption, s=100.00, x=100.00, t=1.00, r=0.10, sigma=0.30, nper_per_year=100, num_trials=1000
   >>> lcall.value()
   27.623315806983594
   >>> lcall.stderr()
   0.7849042509114783
   

   We can obtain a better estimate of the option value (i.e., with lower standard error) by increasing the number of trials:

   >>> # change num_trials to reduce the standard error:
   >>> lcall.num_trials = 100000
   >>> lcall.value()
   28.25988421663273
   >>> lcall.stderr()
   0.07978179212787108
   

   There is no “exact” right answer, but we can increase the number of trials until the standard error is sufficiently small (e.g., maybe a standard error of $0.01 is “good enough”.)

   Here are some additional examples of look-back call options with different parameters:

   >>> lcall = MCLookbackCallOption(35, 30, 1.0, 0.08, 0.25, 100, 100000)
   >>> lcall.value()
   12.77950317693695
   >>> lcall.stderr()
   0.022307140359480636
   >>> lcall = MCLookbackCallOption(35, 40, 1.0, 0.08, 0.40, 100, 100000)
   >>> lcall.value()
   8.794892213218377
   >>> lcall.stderr()
   0.03689429324794948
   

8. Write a class definition for the class MCLookbackPutOption which inherits from the base class MCStockOption. The look-back put option’s payoff the amount by which the option’s exercise price exceeds minimum stock price during the option’s lifetime.

   Calculating the option value:

   The value of the look-back call option (in a single trial) is calculated by:

   where S is price-path of the of the underlying stock form now until time t.

   For example:

   >>> lput = 
   >>> print(lput)
   MCLookbackPutOption, s=100.00, x=100.00, t=1.00, r=0.10, sigma=0.30, nper_per_year=100, num_trials=1000
   >>> lput.value()
   14.968295917890794
   >>> lput.stderr()
   0.3669254592774948
   

   We can obtain a better estimate of the option value (i.e., with lower standard error) by increasing the number of trials:

   >>> # change num_trials to reduce the standard error:
   >>> lput.num_trials = 100000
   >>> lput.value()
   15.173249621749418
   >>> lput.stderr()
   0.03691269486451156
   

   There is no “exact” right answer, but we can increase the number of trials until the standard error is sufficiently small (e.g., maybe a standard error of $0.01 is “good enough”.)

   Here are some additional examples of look-back put options with different parameters:

   >>> lput = MCLookbackPutOption(35, 30, 1.0, 0.08, 0.25, 100, 100000)
   >>> lput.value()
   1.473450162184218
   >>> lput.stderr()
   0.007596162700364556
   >>> lput = MCLookbackPutOption(35, 40, 1.0, 0.08, 0.40,  100, 100000)
   >>> lput.value()
   12.28338712431232
   >>> lput.stderr()
   0.016902847424762495
   

Submitting Your Work

Use the link to GradeScope (left) to submit your work.

Be sure to name your files correctly!

Under the heading for Assignment 7, attach these 3 required files: a7task1.py, a7task2.py and a8graphs.pdf

When you upload the files, the autograder will test your program.

Notes:

• You may resubmit multiple times, but only the last submission will be graded.