Chapter 8

Diffcrcncc cquation for nandom numbcrs. Simple random number generators are based on simulating difference equations. Here is a typical set of two equations: $$ \begin{aligned} &x_{n}=\left(a x_{n-1}+c\right) \bmod m \\ &y_{n}=x_{n} / m \end{aligned} $$ for \(n=1,2, \ldots .\) A seed \(x_{0}\) must be given to start the sequence. The numbers \(y_{1}, y_{2}, \ldots\), represent the random numbers and \(x_{0}, x_{1}, \ldots\) are "help" numbers. Although \(y_{n}\) is completely deterministic from (8.14)(8.15), the sequence \(y_{n}\) appears random. The mathematical expression \(p\) mod \(q\) is coded as \(p \%\) q in Python. Use \(a=8121, c=28411\), and \(m=134456 .\) Solve the system \((8.14)-\) (8.15) in a function that generates and returns \(N\) random numbers. Make a histogram to examine the distribution of the numbers (the \(y_{n}\) numbers are randomly distributed if the histogram is approximately flat). Name of program file: diffeq_random.py.

Independent vs. dependent random numbers. Generate a sequence of \(N\) independent random variables with values 0 or 1 and print out this sequence without space between the numbers (i.e., as 001011010110111010\() .\) The next task is to generate random zeros and ones that are dependent. If the last generated number was 0 , the probability of generating a new 0 is \(p\) and a new 1 is \(1-p\). Conversely, if the last generated was 1, the probability of generating a new 1 is \(p\) and a new 0 is \(1-p\). Since the new value depends on the last one, we say the variables are dependent. Implement this algorithm in a function returning an array of \(N\) zeros and ones. Print out this array in the condense format as described above. Choose \(N=80\) and try the probabilities \(p=0.5, p=0.8\) and \(p=0.9\). Can you by visual inspection of the output characterize the differences between sequences of independent and dependent random variables? Name of program file: dependent_random_variables.py.

Compute \(\pi\) by a random sum. Let \(x_{0}, \ldots, x_{N}\) be \(N+1\) uniformly distributed random numbers between 0 and \(1 .\) Explain why the random sum \(S_{N}=\sum_{i=0}^{N} 2\left(1-x_{i}^{2}\right)^{-1}\) is an approximation to \(\pi\). (Hint: Interpret the sum as Monte Carlo integration and compute the corresponding integral exactly by hand.) Make a program for plotting \(S_{N}\) versus \(N\) for \(N=10^{k}, k=\) \(0,1 / 2,1,3 / 2,2,5 / 2, \ldots, 6\). Write out the difference between \(S_{10^{6}}\) and pi from the math module. Name of program file: MC_pi_plot.py.

\(1 D\) random walk until a point is hit. Set \(\mathrm{np}=1\) in the walk1Dv.py program and modify the program to measure how many steps it takes for one particle to reach a given point \(x=x_{p}\). Give \(x_{p}\) on the command line. Report results for \(x_{p}=\) \(5,50,5000,50000\). Name of program file: walk1Dv_hit_point.py.

Muke u cluss fur \(2 D\) nutuluwr wulk. The purpose of this exercise is to reimplement the ualk2D. py program from Chapter \(8.7 .1\) with the aid of classes. Make a class Particle with the coordinates \((x, y)\) and the time step number of a particle as attributes. A method move moves the particle in one of the four directions and updates the \((x, y)\) coordinates. Another class, Particles, holds a list of Particle objects and a plotstep parameter (as in walk2D.py). A method move moves all the particles one step, a method plot can make a plot of all particles, while a method moves performes a loop over time steps and calls move and plot in each step. Equip the Particle and Particles classes with print functionality such that one can print out all particles in a nice way by saying print \(\mathrm{P}\) (for a Particles instance \(\mathrm{p}\) ) or print self (inside a method). Hint: In _str_-, apply the pformat function from the pprint module to the list of particles, and make sure that__repr__ just reuse __str__ in both classes. To verify the implementation, print the first three positions of four particles in the walk2D. py program and compare with the corresponding results produced by the class-based implementation (the seed of the random number generator must of course be fixed identically in the two programs). You can just perform p.move() and print p three times in a verify function to do this verification task. Organize the complete code as a module such that the classes Particle and Particles can be reused in other programs. The test block should call a run(N) method to run the walk for \(N\) steps, where \(N\) is given on the command line. Compare the efficiency of the class version against the vectorized version in ualk2Dv.py, using the techniques of Appendix G.6.1. Name of program file: walk2Dc.py.

Guess beer brands. You are presented \(n\) glasses of beer, each containing a different brand. You are informed that there are \(m \geq n\) possible brands in total, and the names of all brands are given. For each glass, you can pay \(p\) euros to taste the beer, and if you guess the right brand, you get \(q \geq p\) euros back. Suppose you have done this before and experienced that you typically manage to guess the right brand \(T\) times out of 100 , so that your probability of guessing the right brand is \(b=T / 100\). Make a function simulate \((m, n, p, q, b)\) for simulating the beer tasting process. Let the function return the amount of money earned and how many correct guesses \((\leq n)\) you made. Call simulate a large number of times and compute the average earnings and the probability of getting full score in the case \(m=n=4, p=3, q=6\), and \(b=1 / m\) (i.e., four glasses with four brands, completely random guessing, and a payback of twice as much as the cost). How much more can you earn from this game if your ability to guess the right brand is better, say \(b=1 / 2 ?\) Name of program file: simulate_beer_tasting.py.

Simulate stock prices. A common mathematical model for the evolution of stock prices can be formulated as a difference equation $$ x_{n}=x_{n-1}+\Delta t \mu x_{n-1}+\sigma x_{n-1} \sqrt{\Delta t} r_{n-1} $$ where \(x_{n}\) is the stock price at time \(t_{n}, \Delta t\) is the time interval between two time levels \(\left(\Delta t=t_{n}-t_{n-1}\right), \mu\) is the growth rate of the stock price, \(\sigma\) is the volatility of the stock price, and \(r_{0}, \ldots, r_{n-1}\) are normally distributed random numbers with mean zero and unit standard deviation. An initial stock price \(x_{0}\) must be prescribed together with the input data \(\mu, \sigma\), and \(\Delta t\). We can make a remark that Equation (8.16) is a Forward Euler discretization of a stochastic differential equation for \(x(t)\) $$ \frac{d x}{d t}=\mu x+\sigma N(t) $$ where \(N(t)\) is a so-called white noise random time series signal. Such equations play a central role in modeling of stock prices. Make \(R\) realizations of \((8.16)\) for \(n=0, \ldots, N\) for \(N=5000\) steps over a time period of \(T=180\) days with a step size \(\Delta t=T / N\). Name of program file: stock_prices.py.

Compute with option prices in finance. In this exercise we are going to consider the pricing of so-called Asian options. An Asian option is a financial contract where the owner earns money when certain market conditions are satisfied. The contract is specified by a strike price \(K\) and a maturity time \(T .\) It is written on the average price of the underlying stock, and if this average is bigger than the strike \(K\), the owner of the option will earn the difference. If, on the other hand, the average becomes less, the owner recieves nothing, and the option matures in the value zero. The average is calculated from the last trading price of the stock for each day. From the theory of options in finance, the price of the Asian option will be the expected present value of the payoff. We assume the stock price dynamics given as, $$ S(t+1)=(1+r) S(t)+\sigma S(t) \epsilon(t) $$ where \(r\) is the interest-rate, and \(\sigma\) is the volatility of the stock price. The time \(t\) is supposed to be measured in days, \(t=0,1,2, \ldots\), while \(\epsilon(t)\) are independent identically distributed normal random variables with mean zero and unit standard deviation. To find the option price, we must calculate the expectation $$ p=(1+r)^{-T} \mathrm{E}\left[\max \left(\frac{1}{T} \sum_{t=1}^{T} S(t)-K, 0\right)\right] $$ The price is thus given as the expected discounted payoff. We will use Monte Carlo simulations to estimate the expectation. Typically, \(r\) and \(\sigma\) can be set to \(r=0.0002\) and \(\sigma=0.015 .\) Assume further \(S(0)=100\) a) Make a function that simulates a path of \(S(t)\), that is, the function computes \(S(t)\) for \(t=1, \ldots, T\) for a given \(T\) based on the recursive definition in (8.17). The function should return the path as an array. b) Create a function that finds the average of \(S(t)\) from \(t=1\) to \(t=T\). Make another function that calculates the price of the Asian option based on \(N\) simulated averages. You may choose \(T=100\) days and \(K=102\) c) Plot the price \(p\) as a function of \(N\). You may start with \(N=1000\). d) Plot the error in the price estimation as a function \(N\) (assume that the \(p\) value corresponding to the largest \(N\) value is the "right"price). Try to fit a curve of the form \(c / \sqrt{N}\) for some \(c\) to this error plot. The purpose is to show that the error is reduced as \(1 / \sqrt{N}\). Name of program file: option_price.py. If you wonder where the values for \(r\) and \(\sigma\) come from, you will find the explanation in the following. A reasonable level for the yearly interest- rate is around \(5 \%\), which corresponds to a daily rate \(0.05 / 250=0.0002 .\) The number 250 is chosen because a stock exchange is on average open this amount of days for trading. The value for \(\sigma\) is calculated as the volatility of the stock price, corresponding to the standard deviation of the daily returns of the stock defined as \((S(t+1)-S(t)) / S(t)\). "Normally", the volatility is around \(1.5 \%\) a day. Finally, there are theoretical reasons why we assume that the stock price dynamics is driven by \(r\), meaning that we consider the riskneutral dynamics of the stock price when pricing options. There is an exciting theory explaining the appearance of \(r\) in the dynamics of the stock price. If we want to simulate a stock price dynamics mimicing what we see in the market, \(r\) in Equation \((8.17)\) must be substituted with \(\mu\), the expected return of the stock. Usually, \(\mu\) is higher than \(r \cdot \diamond\)

Compute velocity and acceleration. In a laboratory experiment waves are generated through the impact of a model slide into a wave tank. (The intention of the experiment is to model a future tsunami event in a fjord, generated by loose rocks that fall into the fjord.) At a certain location, the elevation of the surface, denoted by \(\eta\), is measured at discrete points in time using an ultra-sound wave gauge. The result is a time series of vertical positions of the water surface elevations in meter: \(\eta\left(t_{0}\right), \eta\left(t_{1}\right), \eta\left(t_{2}\right), \ldots, \eta\left(t_{n}\right)\). There are 300 observations per second, meaning that the time difference between to neighboring measurement values \(\eta\left(t_{i}\right)\) and \(\eta\left(t_{i+1}\right)\) is \(h=1 / 300\) second. Write a Python program that accomplishes the following tasks: 1\. Read \(h\) from the command line. 2\. Read the \(\eta\) values in the file src/random/gauge. dat into an array eta. 3\. Plot eta versus the time values. 4\. Compute the velocity \(v\) of the surface by the formula $$ v_{i} \approx \frac{\eta_{i+1}-\eta_{i-1}}{2 h}, \quad i=1, \ldots, n-1 $$ Plot \(v\) versus time values in a separate plot. 5\. Compute the acceleration \(a\) of the surface by the formula $$ a_{i} \approx \frac{\eta_{i+1}-2 \eta_{i}+\eta_{i-1}}{h^{2}}, \quad i=1, \ldots, n-1 $$ Plot \(a\) versus the time values in a separate plot.

Numerical differentiation of noisy signals. The purpose of this exercise is to look into numerical differentiation of time series signals that contain measurement errors. This insight might be helpful when analyzing the noise in real data from a laboratory experiment in Exercises \(8.44\) and \(8.46\). 1\. Compute a signal $$ \bar{\eta}_{i}=A \sin \left(\frac{2 \pi}{T} t_{i}\right), \quad t_{i}=i \frac{T}{40}, i=0, \ldots, 200 $$ Display \(\bar{\eta}_{i}\) versus time \(t_{i}\) in a plot. Choose \(A=1\) and \(T=2 \pi\). Store the \(\bar{\eta}\) values in an array etabar. 2\. Compute a signal with random noise \(E_{i}\), $$ \eta_{i}=\bar{\eta}_{i}+E_{i} $$ \(E_{i}\) is drawn from the normal distribution with mean zero and standard deviation \(\sigma=0.04 A\). Plot this \(\eta_{i}\) signal as circles in the same plot as \(\eta_{i}\). Store the \(E_{i}\) in an array E for later use. 3\. Compute the first derivative of \(\bar{\eta}_{i}\) by the formula $$ \frac{\bar{\eta}_{i+1}-\bar{\eta}_{i-1}}{2 h}, \quad i=1, \ldots, n-1 $$ and store the values in an array detabar. Display the graph. 4\. Compute the first derivative of the error term by the formula, $$ \frac{E_{i+1}-E_{i-1}}{2 h}, \quad i=1, \ldots, n-1 $$ and store the values in an array \(d E .\) Calculate the mean and the standard deviation of \(\mathrm{dE}\). 5\. Plot detabar and detabar \(+\mathrm{dE}\). Use the result of the standard deviation calculations to explain the qualitative features of the graphs. 6\. The second derivative of a time signal \(\eta_{i}\) can be computed by $$ \frac{\eta_{i+1}-2 \eta_{i}+\eta i-1}{h^{2}}, \quad i=1, \ldots, n-1 $$ Use this formula on the etabar data and save the result in d2etabar. Also apply the formula to the E data and save the result in \(\mathrm{d} 2 \mathrm{E}\). Plot d2etabar and d2etabar + d2E. Compute the standard deviation of d2E and compare with the standard deviation of \(\mathrm{dE}\) and \(\mathrm{E} .\) Discuss the plot in light of these standard deviations. Name of program file: sine_noise.py. 0