SciPy library ODE solver

1    SciPy library ODE solver (20% of project credit)

 

1.1    Overview

The recommended SciPy ODE solver is solve ipv. If one needed to numerically solve an ODE in Python, this would be a standard library function to use. One would be expected to read the documentation as necessary in order to use the function correctly. This part of the project is designed to let you do just this.

The Lorenz model is given by the following system of ODEs

 

= β(3 . y)]

dt

d3

d之

= y3 . ←之]

dt

(1)

 

(2)

 

(3)

where β , β and ← are parameters. The goad of this part of the project is to use solve ipv to solve this system of ODEs according to the specifications below.

solve ipv is similar to, but also distinctly different from odeint, and you will need to understand these differ- ences and account for them.

The Wikipedia page on the Lorenz system provides useful overview together with many plots.

 

1.2    Particulars

● The parameters values are: β = 10, β = 28 and ← = 8|3.

● The initial condition (y(0)] 3(0)] 之(0)) is your choice. Explore different possibilities and then choose one.

● The time interval is [0] tf ], where the final time tf  is your choice. You will want tf  larger than 10 but smaller than 100.

● In all cases you will need the input parameterfun, t span and y0. Use default values of input parame- ters unless required to do otherwise.

 

1.3    Computational Tasks

● After reading the SciPy documentation, write a Python function that computes and returns the right- hand side (RHS) of the Lorenz system of ODEs in a form required by solve ipv.

● Write Python code to call solve ipv, and then extract and plot the solution. You should be able to plot time series: y(t), 3(t), 之 (t), as function of t, and phase portraits: standard ones would be (y] 3) and (y] 之).  You may produce 3D phase portraits if you wish.  For the time series you can plot all three variables on one plot or on multiple plots as you think best.

The module did not cover the return type of solve ipv, but it is not difficult to obtain the needed solution from the return from solve ipv.

● You will find that to get smooth plots, the input argument t eval must be used. Understand it and set values that give you smooth looking plots.

● Run cases for different values of the initial condition and final time tf  to decide values that you will use for the initial condition and final time. Then leave this fixed for the remaining tasks.

● Run and plot the solution using the time points returned by solve ipv when t eval is not included as an input parameter.

● Run and plot the solution using t eval as an input parameter set with reasonable values to give smooth plots.

● Run the following cases and report that values (y] 3] 之) at the final time. You do not need to produce plots for these cases:

–  Default parameters.

–  Default parameters, except the method is explicit Runge-Kutta method of order 8, rather than explicit Runge-Kutta method of order 5(4).

–  Default parameters, except the relative tolerance of the solver is 10 −6  rather than the default 10 −3 .

–  Default parameters, except the method is explicit Runge-Kutta method of order 8 and the relative tolerance of the solver is 10 −6 .

 

1.4    Report contents

See general discussion of report contents in Sec. 5 and 6.

 

● A short introduction should state the equations begin solved, the parameter being used, the initial condition and time interval for your solutions and plots.

● A section computing and showing time series and phase portraits, first with using the time points returned by solve ipv and then again using t eval.  The exact content of the plots is your choice, including whether to show any 3D plots. Make them nice.

● A section computing and reporting final values for the four cases run with different methods and relative tolerances. Are the solutions converged? Look up “sensitive dependence on initial conditions” and write a short comment (a few sentences) about this. What are the implications for examining the accuracy of your ODE solution?

 

For reference, with the code cells collapsed, the part of the project should be roughly the equivalent of 2 sides of A4, depending on the size and number of your plots.

2    Pricing a Basket Option (40% of project credit)

 

2.1    Overview

A basket option is a financial derivative whose underlying is a weighted sum or average of different assets sm (t) that have been grouped together in a basket. We will only consider the case where the basket is the average of the sm (t).

M

l(t) =         sm (t)                                                           (4)

 

Options differ in their payoff at expiry t = r. We will consider

 

 

● European call option whose payoff is

max(l(r) . D] 0)]

where D is the strike price.

● Asian call option whose payoff is

max( . D] 0)]

 

(5)

 

 

 

(6)

where  is the mean value of l(t) over the life of the option 0 < t < r, and D is the strike price.

 

 

● Lookback call option whose payoff is max(l(r) . b] 0)] where b is the minimum value of l(t) over the life of the option 0 < t < r. In each case, the option price is the expected discounted payoff. We assume the underlying processes for the sm are geometric Brownian motion dsm (t) = rsm (t)dt + βm sm (t)dwm (t)]         where r and βm are constants and where the Wiener processes wm are correlated E[dwi (t)dwj (t)] = βij dt   or    E[βi dwi (t)βj dwj (t)] = Σij dt] where βij is the correlation matrix and Σij is the covariance matrix, assumed to be constant in time. The aim of this part of the project is to price basket options by Monte Carlo methods.

(7)               (8)

 

2.2    Particulars

● The basket is composed of k = 5 assets: s1 (t)] · · · ] s5 (t).

● l(t) is given by Eq. (4), the average of the sm (t).

● The correlation matrix β is given by:

 

 

┌' 01σ1 β =  '(') .0 σ 1 '  0   ● The volatilities βk are given by:

0 σ 1 1 0 0 0 σ 2

.0 σ 1 0   1   0   0

0 0 0 1 0 σ 15

00σ2 ┐' 0   '(') 0 σ 15 '

 

 

β 1  = 0 σ 1] β2  = 0 σ 12] β3  = 0 σ 13] β4  = 0 σ 09] β5  = 0 σ 11

● The spot prices of the assets are:

s1 (0) = 100] s2 (0) = 90] s3 (0) = 85] s4 (0) = 105] s5 (0) = 120

● The interest rate is r = 0 σ 05.

● Time r to maturity is 3 years.

● The strike price is D = 100, which make it at-the-money since l(0) = 100.

● The mean price  and minimum price b are based on daily closings, so

 

N

 =        l(tn )

 

b = min l(tn )

where tn  = ndt. The time step dt is one day where we assume there are 260 (trading) days in a year.

● Do not exceed Npaths  = 104 Monte Carlo paths in your computations!

 

2.3    Computational Tasks

● Implement a Python function that performs Monte Carlo simulations of Eqs. (8) using Euler time stepping with a time step of one day.  The function should have input parameters: a vector of spot prices s0, as well as r, r, Σ , the number of paths to compute Npaths , and any other necessary parameters. The function should return the time grid t and a three-dimensional array smnp  of daily closings, where m is the asset label (1 to 5), n is the time step, and p the path index.

● Write Python code that sets up the problem parameters, calls the Monte Carlo function, and processes the returned array smnp  to obtain the three option prices, together with 95% confidence intervals based on the number of Monte Carlo paths, Npaths , used in your simulations.

–  It is suggested that you first consider only the European call option. Test your code by setting Σ to a diagonal matrix. This effectively reduces the computations to the cases of independent European call options. Test your Monte Carlo code against the Black Scholes formula until you are satisfied that it is correct. You may leave in testing code in your submitted project, but you should comment it out.

–  Continue by implementing the other payoff functions and reinstate the covariance matrix and initial asset prices as specified above.

● Implement a Python function that performs Monte Carlo simulations of Eqs. (8) using antithetic vari- ance reduction. All input parameters should be the same as the previous Monte Carlo function. Put this function in the same code cell as the previous function.

● Write Python code that calls the antithetic Monte Carlo function and computes option prices, together with 95% confidence intervals, for the three option types.

 

2.4    Report contents

See general discussion of report contents in Sec. 5 and 6.

 

● The focus of this part of the project is Python coding. Your Python code should be clearly structured and contain sufficient comments so that the purpose of each part of the code is clear to a Python programmer. Use the structure of the Jupyter notebook to block your code into clearly defined units.

● The report should address the improvements of the variance reduction.  There are three good ap- proaches here.

The first is a raw comparison of the 95% confidence intervals for a fixed sample size Npaths with and without variance reduction.

The second is to compare the number of paths that would be needed to obtain a specific absolute accuracy with and without variance reduction.

The third is a comparison of histograms of the discounted payoffs with and without variance reduction.

● The report should discuss the problem with using antithetic variance reduction for the lookback option. Include plots if you find that helpful.

● If you were able to resolve the problem and price the lookback option, explain how.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6

3    Machine Learning:  Bike Sharing Demand Data (20% of project credit)

 

3.1    Overview

In a recent research article published in the journal Computer Communications1 , Sathishkumar et al. sought to predict the ”bike count required at each hour for the stable supply of rental bikes”.  They employed a number of regression models. The goal of this part of the project is to reproduce some aspects of this study. You will have the opportunity to compare your machine-learning results with a published study and you will have the opportunity to use pandas and seaborn to generate various plots.

 

3.2    Particulars

A link to the original research article and a modified dataset are on the Module Resources web page where weekly Jupyter notebooks appear. You will need to refer to the article for some of the tasks below. You will need to download SeoulBikeData mod.csv and put it into the folder with your project notebook.

SeoulBikeData mod.csv has been modified from the original dataset to remove the categorical variables, and to convert dates to months. Months have been coded by number, e.g.  1 = January, etc.  Only the first six months are included in the modified dataset.

Scikit-learn will be used for all machine learning tasks. Pandas and seaborn will be used for inspecting the data and making plots.

 

3.3    Tasks

● Using pandas, read SeoulBikeData mod.csv into a Dataframe and ‘describe‘ the Dataframe.

● Plot a histogram of ‘Rented Bike Count‘ and make a box plot similar that shown in Fig. 3 of Ref. [1]. Try to match the aspect ratios of the published figures.

● Produce two violin plots: one showing ‘Rented Bike Count‘ for different values of the ‘Month‘ and the other showing ‘Rented Bike Count‘ for different values of the ‘Hour‘.

● From the full Dataframe, create a design matrix x (as a Dataframe) containing all the columns except ‘Rented Bike Count‘ and the target 3 (as a Series) containing only the ‘Rented Bike Count‘ column.

● Perform a test-train split. You must use the same percentage of data for testing and training as was used in Ref. [1], and you must state what percentages you are