So far all of our methods have assumed that the abscissa points have been equally spaced along the $x$-axis

We could have changed this if we wanted but there seemed little point and it would have made the algebra more complicated

The question therefore arises: Can we improve the speed or accuracy of our routine by using a different set of abscissa points?

The answer is YES and the set of methods we can use are referred to as Gaussian quadrature

In general all quadrature methods are based on the formula:

$\displaystyle\int_{a}^{b}f(x)dx\approx \displaystyle\sum_{i=1}^{n}\omega_i f(x_i)$

Where $\omega_i$ are a set of weights and $x_i$ are a set of abscissa points

As we are now free to choose $2n$ values we should be able to approximate our integral much more accurately - in fact if we motivate our selection of $\omega_i$ and $x_i$ by fitting polynomials then this method will calculate the integral of an order $2n-1$ polynomial exactly

W.l.o.g we work on the interval $[-1, 1]$

We wish to find $\omega_1$, $\omega_2$, $x_1$ and $x_2$ such that

$\displaystyle\int_{-1}^{1}f(x)dx \equiv \displaystyle\sum_{i=1}^{2}\omega_i f(x_i)$, for a cubic $f(x)$, i.e.:

$\displaystyle\int_{-1}^{1} a_0+a_1x+a_2x^2+x_3x^3dx \equiv \omega_1 f(x_1) + \omega_2 f(x_2)$, i.e. this relation must hold for all choices of $a_i$

We proceed with basic calculus

$\left[ a_0x+\frac{a_1x^2}{2}+\frac{a_2x^3}{3}+\frac{x_3x^4}{4} \right]_{-1}^{1} \equiv \omega_1 f(x_1) + \omega_2 f(x_2)$

$a_0(1-(-1)) + a_1\frac{1^2-(-1)^2}{2} + a_2\frac{1^3-(-1)^3}{3} + a_3\frac{1^4-(-1)^4}{4} \equiv \omega_1 f(x_1) + \omega_2 f(x_2)$

$2a_0 + \frac{2}{3}a_2 \equiv \omega_1 \left(a_0+a_1x_1+a_2x_1^2+a_3x_1^3\right)+ \omega_2\left(a_0+a_1x_2+a_2x_2^2+a_3x_2^3\right)$

$2a_0 + \frac{2}{3}a_2 \equiv a_0(\omega_1 + \omega_2) + a_1(\omega_1 x_1 + \omega_2 x_2) + a_2(\omega_1 x_1^2 + \omega_2 x_2^2) + a_3(\omega_1 x_1^3 + \omega_2 x_2^3)$

As this is an identity we can match the co-efficients of the $a_i$

$a_0$: $2 = \omega_1 + \omega_2$

$a_1$: $0 = \omega_1 x_1 + \omega_2 x_2$

$a_2$: $\frac{2}{3} = \omega_1 x_1^2 + \omega_2 x_2^2$

$a_3$: $0 = \omega_1 x_1^3 + \omega_2 x_2^3$

The best way to solve these equations is guess and check as follows:

as symmetry seems likely assume $\omega_1 = \omega_2$ therefore $\omega_1 = \omega_2 = 1$

$a_1$ gives $x_1 = -x_2$

$a_2$ gives $\frac{2}{3} = \omega_1 x_1^2 + \omega_2 x_2^2 = 2 x_1^2$

so $3x_1^2-1=0$ - this is the Legendre Polynomial order 2

w.l.o.g. $x_1=-\sqrt{\frac{1}{3}}$ and $x_2=\sqrt{\frac{1}{3}}$

Calculate $\displaystyle\int_{-1}^{1} 2+6x-3x^2+x^3dx$, both analytically and using 2 point Gaussian quadrature

The last thing we need to do to make our method useful is to drop the requirement to work on the $[-1, 1]$ interval

There are two additional things we need to consider to to this:

- the $\omega_i$ need to be transposed from the [-1,1] interval to the [5,10] interval
- The area calculated at the end needs to be ratioed up from a width of 2 (1 - (-1)) to a width of 5 (10-5)

Calculate $\displaystyle\int_{5}^{10} 2+6x-3x^2+x^3dx$, both analytically and using 2 point Gaussian quadrature

We could repeat the above exercise for 3 and then 4,5 ,6 etc points. The algebra does get complicated, but fortunately we have the following theorem:

If we perform n point quadrature by selecting the $x_i$ to be the roots of the Legendre Polynomials and the weighting function to be $\omega_i = \frac{2}{(1-x_i^2)(P'_n(x_i))^2}$, where $P'_n(x_i)$ is the differential of the $n^{th}$ Legendre Polynomial then our resulting quadrature will integrate polynomials of order $2n-1$ exactly

We do not attempt to prove this result here - but the above analysis is the proof of this result for 2 point quadrature. You can see how the Legendre Polynomial arises in the above analysis

The Legendre Polynomials are defined by: $P_n(x)=\frac{1}{2^n} \displaystyle\sum_{k=0}^{n} \left({n \choose k}^2 \times (x-1)^{n-k} \times (x+1)^k\right)$

Calculating $P'_n(x)$ is a simple application of the product rule

So we can see that the first few of the Legendre Polynomials are:

- $n=0: P_0(x)=0$
- $n=1: P_1(x)=x$
- $n=2: P_2(x)=\frac{1}{2} \left(3x^2-1 \right)$
- $n=3: P_3(x)=\frac{1}{2} \left(5x^3-3x \right)$
- $n=4: P_4(x)=\frac{1}{8} \left(35x^4-30x^2 + 3 \right)$

Produce a spreadsheet which draws graphs of each of the Legendre Polynomials

What do you notice about the dispersion of the roots?

Now compare the roots of the n-th order Legendre Polynomial with the value of $x=cos \left( \frac{(i-0.25) \pi}{n+0.5} \right)$

How might we use the formula: $x=cos \left( \frac{(i-0.25) \pi}{n+0.5} \right)$

To produce initial estimates of the roots so we can then use Newton Rhapson to produce our exact roots

1. Develop a VBA function which can perform an n point Gaussian quadrature over a given interval for any given function

2. change your code so that the calculation of the Legendre Polynomial roots is only calculated once, however many times the routine is called

Hint: you will need to store the $\omega_i$ and $x_i$ values outside of the function so that they have module wide scope and use a boolean array to determine whether the values have already been calculated or not

3. Use 2, 4, 6, 8 and 10 point Gaussian quadrature to calculate $\displaystyle\int_{0}^{\pi} sin x dx$

To use complex numbers in VBA we make use of the OOP (object orientated programming) features

VBA is not a natural OOP language is the same way that Delphi or C++ are, but objects can still be used

To see how OOP works we are going to look at manipulating complex numbers. Enter the following code into your Module

Type cnum x As Double y As Double End Type Function add_cnum(ByRef one As cnum, ByRef two As cnum) As cnum add_cnum.x = one.x + two.x add_cnum.y = one.y + two.y End Function Function set_cnum(ax, ay As Double) As cnum set_cnum.x = ax set_cnum.y = ay End Function Sub do_calc() Dim c1 As cnum Dim c2 As cnum Dim c3 As cnum c1 = set_cnum(Cells(1, 1), Cells(1, 2)) c2 = set_cnum(Cells(2, 1), Cells(2, 2)) c3 = add_cnum(c1, c2) Cells(3, 1) = c3.x Cells(3, 2) = c3.y End Sub

Enter two complex numbers into the cells A1:B1 and A2:B2

Click the Run Macro button in Excel, select *do_calc* and you will see the two complex numbers added together in the cells A3:B3

There are a number of features of this code:

- A new type (or object) cnum is created
- This new type has two variables within it: a and b and so is in total a complex number
- This new 'object' can be passed as a parameter to functions and returned as a result from functions.

Write new functions: mult_cnum, div_cnum, exp_cnum, ln_cnum and power_cnum

The following algebra will help:

$(x_1 + iy_1) \times (x_2 + iy_2) = x_1 x_2 - y_1 y_2 + i(y_1 x_2 + x_1 y_2)$

To divide two complex numbers you will need to use the complex conjugate

$\frac{1}{x+iy} = \frac{x-iy}{(x+iy)\times(x-iy)} = \frac{x-iy}{x^2+y^2}$

To calculate the exponential you will need Euler's Formula:

$e^{ix}=cos x + i sin x$

To calculate the logarithm you will need to make $x_2$ and $y_2$ the subjects of the formula, starting from: $x_1 + iy_1 = e^{x_2 + iy_2}$ and again applying Euler's formula and matching terms in real and imaginary parts

To calculate a power you need to start with $a=b^c$ and follow the logic below

$$a=b^c$$

$$ln(a) = c\times ln(b)$$

$$a=exp(c\times ln(b))$$

Fourier analysis allows us to convert waveforms into spectral densities and back again

In excel draw a graph of $y=sin x$

On the same axes draw $y=sinx + \frac{sin 2x}{2}$

Create a VBA function for $y=\displaystyle\sum_{k=1}^{n} \frac{sin kx}{k}$ and again plot on the same axes for different values of $n$

You can soon see that as $n \rightarrow \infty$ this function becomes a sawtooth wave

This raises the question: If we started with the sawtooth wave wave could we calculate the co-efficients

This is where fourier analysis comes in

For our purposes we will skip Fourier series and go onto the continuous Fourier transform

The fourier tranform $\hat{f}$ of an integrable function $f:\Re \rightarrow C$ is given by

$\hat{f}(\xi)=\displaystyle\int_{-\infty}^{\infty}f(x) e^{-2 \pi ix \xi} dx$

$f$ is the original wave function of time $x$ and $\hat{f}$ is the spectral density of frequency $\xi$

The process can be inverted by the inverse transform:

$f(x)=\displaystyle\int_{-\infty}^{\infty} \hat{f}(\xi) e^{2 \pi ix \xi} d\xi$

For our purposes (later in the course when we look at stochastic volatility) we use Fourier inversion to get from the characteristic function of a probability distribution back to the probability distribution

The above formulas look rather daunting but we start by breaking them down with Euler's formula:

$e^{ix}=cos x + i sin x$

Write a VBA function which takes a function $f$, a parameter $\xi$, specifies bounds $a$ and $b$ outside of which the function will be close to zero, and returns the Fourier transform $\hat{f}(\xi)$

Write another VBA function which takes a function $\hat{f}$, a parameter $x$, specifies bounds $a$ and $b$ outside of which the function will be close to zero, and returns the inverse Fourier transform $f(x)$

Test your VBA code of the function: $f(x)=e^{-x^2}$

i.e. check if the Fourier transform followed by the inverse transform returns you to the original function

Investigating whether or not the value calculated above actually makes any sense will form part of your written coursework. For the calculation part of the coursework it is sufficient to be able to just calculate what the actual number is.

Many complicated operations can be simplified with matrix notation. A suit of VBA routines which can handle basic matrix manipulation will make many tasks much easier

You have already written a routine that multiplies two square matrices together

$$ \left( \begin{array}{ccc} a_{11} & a_{12} & ... \\ a_{21} & ... & ... \\ ... & ... & ... \end{array} \right) \times \left( \begin{array}{ccc} b_{11} & b_{12} & ... \\ b_{21} & ... & ... \\ ... & ... & ... \end{array} \right) = \left( \begin{array}{ccc} a_{11}b_{11} + a_{12}b_{21} + ... & ... & ... \\ a_{21}b_{11} + a_{22}b_{21}+... & ... & ... \\ ... & ... & ... \end{array} \right)$$

Write a function which can take two general matrices and multiply them together returning the answer to the calling routine.

Write another sub behind an excel button which calls your function and multiplies two excel matrices together putting the output back on the excel sheet. (This is to test your function)

Before working out the inverse of a matrix we need to calculate the determinant. The determinant of a 2 by 2 matrix is simple and you can see how the pattern develops for higher order matrices

$${\begin{vmatrix}a&b\\c&d\end{vmatrix}} = \frac{1}{ad-bc}$$

To calculate the determinant of a 3 by 3 matrix we have to use the determinants of the 2 by 2 matrix

$${\begin{vmatrix}a&b&c\\d&e&f\\g&h&i\end{vmatrix}} =a{\begin{vmatrix}e&f\\h&i\end{vmatrix}}-b{\begin{vmatrix}d&f\\g&i\end{vmatrix}}+c{\begin{vmatrix}d&e\\g&h\end{vmatrix}}$$

This pattern then continues for matrices of any order

$${\begin{vmatrix}a&b&c&d\\e&f&g&h\\i&j&k&l\\m&n&o&p\end{vmatrix}} =a{\begin{vmatrix}f&g&h\\j&k&l\\n&o&p\end{vmatrix}} -b{\begin{vmatrix}e&g&h\\i&k&l\\m&o&p\end{vmatrix}}+c{\begin{vmatrix}e&f&h\\i&j&l\\m&n&p\end{vmatrix}}-d{\begin{vmatrix}e&f&g\\i&j&k\\m&n&o\end{vmatrix}}$$

Write a function which calculates the determinant of a $n \times n$ matrix

being able to invert a matrix is very useful as it enables us to directly solve systems of linear equations

The method is as follows:

- Calculate the cofactor matrix
- Transpose the cofactor matrix to get the Adjoint matrix
- Divide the adjoint matrix by the determinant to get the inverse matrix

The cofactor matrix is formed by taking each element and then taking the determinant of the matrix which is left when you remove the row and column of that element from the matrix. This value is then negated if the sum of the row and column is odd - s0 that the signs alternate

So if we take the matrix $\left( \begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{array} \right)$ then the coFactor of $a_{21}$ is minus the determinant of $\left( \begin{array}{cc} a_{12} & a_{13} \\ a_{32} & a_{33} \end{array} \right)$

Given the matrix: $\left( \begin{array}{ccc} 4 & 7 & 1 \\ 3 & 3 & 8 \\ 12 & 1 & 11 \end{array} \right)$, Calculate its inverse

The matrix of determinants: $\left( \begin{array}{ccc} 25 & 63 & -33 \\ -76 & 32 & 80 \\ 53 & -29 & -9 \end{array} \right)$

The is the transpose of the cofactor matrix: $\left( \begin{array}{ccc} 25 & -76 & 53 \\ 63 & 32 & -29 \\ -33 & 80 & -9 \end{array} \right)$

${\begin{vmatrix}4 & 7 & 1\\ 3 & 3 & 8 \\ 12 & 1 & 11 \end{vmatrix}} =4{\begin{vmatrix} 3 & 8 \\ 1 & 11 \end{vmatrix}}-7{\begin{vmatrix} 3 & 8 \\ 12 & 11 \end{vmatrix}}+1{\begin{vmatrix} 3 & 3 \\ 12 & 1 \end{vmatrix}}$

$=4 \times 25 - 7 \times (-63) + 1 \times (-33) = 508$

$A^{-1}=\frac{1}{508} \times \left( \begin{array}{ccc} 25 & -76 & 53 \\ 63 & 32 & -29 \\ -33 & 80 & -9 \end{array} \right)$

How can you check if your inverse matrix is correct

Check that $A \times A^{-1} = I$

If you wish to fit a function of several parameters to a set of data - this can be a very computationally intensive task

If there is only one parameter then using goal seek in Excel will probably be perfectly serviceable, but you cannot use this method if you have multiple parameters

first you need to decide what is meant by best fit. For most purposes minimising the square of the deviation between the actual data and the fitted data will be a sensible metric, hence the expression "Least squares Regression"

Then we need to decide how to search the vector space of all possible parameter combinations for the one for which the sum of the least square deviation is the least

We list below the a number of methods in order of increasing sophistication. We assume:

- there are $n$ parameters
- the sum least squares function is denoted by $f$
- the current set of parameters is stored in a vector denoted by $P$

The methods are:

- Calculate $f$ at every point in a $n$ dimensional grid and then select the point in the grid for which $f$ is least. This is the set of parameters $P$ which gives the best fit
- Minimise for parameter 1, keeping the other parameters constant, then do the same for parameter 2 and so on, starting again with parameter 1 and repeating until the function cannot be reduced by changing any further parameters
- Calculate the gradient vector at $P_0$ and then minimise in the steepest direction. Once the minimum has been found in this direction recalculate the gradient and repeat the process until you reach a minimum in all directions

On first inspection the second and third methods would appear to be quite good methods, but in practice they both tend to zig zag very slowly towards a solution and a much more efficient method is found by the method of conjugate gradients

Before we set out the conjugate gradient methods we need to understand how to do an individual line minimisation as which ever method we use (apart from the grid) we will need to minimise out function along any particular line in the vector space of the parameters

Line minimisation is analogous to Newton-Rhapson root finding although we are not solving $f(x)=0$ but rather we are solving $f'(x)=0$ i.e. the point which the function (least squares deviation) is at a minimum

For Newton Rhapson the assumption we made was of local linearity - sop this time the assumption we make is that the function is locally quadratic

We can use the first and second derivatives if we have them (a vector and matrix respectively) but if not a simple numerical procedure will suffice

calculate $f_0=f(\textbf{P}), f_1+f(\textbf{P}+h\textbf{x}), f_2=f(\textbf{P}+2h\textbf{x})$, where $\textbf{x}$ is the unit vector in parameter space of the direction in which we are minimising and $h$ is a small increment

Then $f''_0 = \frac{(f_2-f_1)-(f_1-f_0)}{h^2}$ is the rate at which the (negative) gradient is getting shallower and so given $f'_0 = \frac{(f_1-f_0)}{h}$ is the current gradient in the direction $\textbf{x}$ then our initial line minimisation guess is $P_{n+1} = P_n - \textbf{x} \times \frac{f'_0}{f''_0}$

We then repeat this process along the same line until we have found our one dimensional minimum - at which point we continue with whatever multi-dimensional method we were using

The conjugate gradient method is as follows:

Step 1: Initialise the set of directions $\textbf{u}_i$ to the basis vectors $\textbf{e}_i$

Step 2: Save your starting position as $\textbf{P}_0$

Step 3: For $i=1,...,N$, move $\textbf{P}_{i-1}$ to the minimum along direction $\textbf{u}_i$ and call this point $\textbf{P}_i$

Step 4: For $i=1,...,N$, set $\textbf{u}_i$ to be $\textbf{u}_{i+1}$

Step 5: Set $\textbf{u}_N$ to be $\textbf{P}_N - \textbf{P}_0$

Step 6: Move $\textbf{P}_N$ to the minimum along the direction $\textbf{u}_N$ and call this point $\textbf{P}_0$

Step 7: Repeat this basic algorithm $N$ times and any quadratic form will be minimised

The result in step 7 was proved by Powell in 1964

This spreadsheet Powell.xls contains an implementation of the method

Generating functions are a strange concept when you first come across them but they have many applications and often make otherwise complicated calculations much more straighforward

The probability generating function of a random variable $X$ is $G(t)=E \left(t^X \right)$

For a discrete random variable (on the integers) this is: $G(t)=\sum\limits_{x=-\infty}^{\infty} p(x) t^x$

For a continuous random variable this is: $G(t)=\int\limits_{-\infty}^{\infty} f(x) t^x dx$, where $f(x)$ is the probability density function.

What is the probability generating function which describes a dice roll?

$G(t)=\frac{1}{6}\left(t^1 + t^2 + t^3 + t^4 + t^5 + t^6 \right)$

Now calculate $\left(G(t)\right)^2$

Can you see the practical applciation that probability generating function might have?

You can read off the co-efficients and they give the probability of each score when you add two dice rolls together

The moment generating function of a random variable $X$ is $M(t)=E \left(e^{tX} \right)$

For a continuous random variable this is: $M(t)=\int\limits_{-\infty}^{\infty} f(x) e^{tx} dx$, where $f(x)$ is the probability density function.

Calculate the moment generating function for the the normal distribution:

for $X\sim N(\mu, \sigma^2), M_X(t)=e^{t\mu+\frac{1}{2}\sigma^2t^2}$

Now differentiate $M_X(t)$ and find $\frac{dM_X(0)}{dt}$ and $\frac{d^2M_X(0)}{dt^2}$

$\frac{dM_X(0)}{dt} = \mu$

$\frac{d^2M_X(0)}{dt^2} = \mu + \sigma^2$

The cumulant generating function of a random variable $X$ is $K(t)=ln \left( E \left(e^{tX} \right)\right)$

We can see that for $X\sim N(\mu, \sigma^2), K_X(t)=t\mu+\frac{1}{2}\sigma^2t^2$

So $K'(t)=\mu + \sigma^2 t$ and $K''(t)=\sigma^2$

The cumulants are defined by $k_n=\frac{d^nK(0)}{dt^n}$

The cumulants are *$k_1=\mu$, $k_2=\sigma^2$, $k_n=0$ for $n \geqslant 3$*