In this chapter we introduce some of the mathematical concepts that will be needed to deal with the option pricing and stochastic volatility models introduced in this book, and to help readers implement these concepts as functions and routines in VBA. First, we introduce complex numbers, which are needed to evaluate characteristic functions of distributions driving option prices. These are required to evaluate the option pricing models of Heston (1993) and Heston and Nandi (2000) covered in Chapters 5 and 6, respectively. Next, we review and implement Newton’s method and the bisection method, two popular and simple algorithms for finding zeros of functions. These methods are needed to find volatility implied from option prices, which we introduce in Chapter 4 and deal with in Chapter 10. We show how to implement multiple linear regression with ordinary least squares (OLS) and weighted least squares (WLS) in VBA. These methods are needed to obtain the deterministic volatility functions of Chapter 4. Next, we show how to find maximum likelihood estimators, which are needed to estimate the parameters that are used in option pricing models. We also implement the Nelder-Mead algorithm, which is used to find the minimum values of multivariate functions and which will be used throughout this book. Finally, we implement cubic splines in VBA. Cubic splines will be used to obtain model-free implied volatility in Chapter 11, and model-free skewness and kurtosis in Chapter 12.
Most of the numbers we are used to dealing with in our everyday lives are real numbers, which are defined as any number lying on the real line ℜ = (−∞, + ∞). As such, real numbers can be positive or negative; rational, meaning that they can be expressed as a fraction; or irrational, meaning that they cannot be expressed as a fraction. Some examples of real numbers are 
and π. Complex numbers, however, are constructed around the imaginary unit i defined as 
. While i is not a real number, i2 is a real number since i2 = −1. A complex number is defined as a = x + iy, where x and y are both real numbers, called the real and imaginary parts of a, respectively. The notation Re[] and Im[] is used to denote these quantities, so that Re[a] = x and Im[a] = y.
Many of the operations on complex numbers are done by isolating the real and imaginary parts. Other operations require simple tricks, such as rewriting the complex number in a different form or using its complex conjugate. Krantz (1999) is a good reference for this section.
Addition and subtraction of complex numbers is performed by separate operation on the real and imaginary parts. It requires adding and subtracting, respectively, the real and imaginary parts of the two complex numbers:
Multiplying two complex numbers is done by applying the distributive axiom to the product, and regrouping the real and imaginary parts:
The complex conjugate of a complex number is defined as 
and is useful for dividing complex numbers. Since 
, we can express division of any two complex numbers as the ratio
Exponentiation of a complex number is done by applying Euler’s formula, which produces
Hence, the real part of the resulting complex number is exp(x) cos(y), and the imaginary part is exp(x) sin (y). Obtaining the logarithm of a complex number requires algebra. Suppose that w = a + ib and that its logarithm is the complex number z = x + iy, so that z = log(w). Since w = exp(z), we know that a = ex cos (y) and b = ex sin (y). Squaring these numbers, applying the identity cos (y)2 + sin (y)2 = 1, and solving for x produces 
. Taking their ratio produces
and solving for y produces y = Im[z] = arctan(b/a).
It is now easy to obtain the square root of the complex number w = a + ib, using DeMoivre’s Theorem:
By arguments in the previous paragraph, we can write w = r cos (y) + ir sin (y) = reiy, where y = arctan(b/a) and 
The square root of w is therefore
Applying DeMoivre’s Theorem with n = 1/2, this becomes
so that the real and imaginary parts of 
, respectively.
Finally, other functions of complex numbers are available, but we have not included VBA code for these functions. For example, the cosine of a complex number z = x + iy produces another complex number, with real and imaginary parts given by cos (x) cosh (y) and − sin (x) sinh (y) respectively, while the sine of a complex number has real and imaginary parts sin (x) cosh (y) and − cos (x) sinh (y), respectively. The hyperbolic functions cosh (y) and sinh (y) are defined in Exercise 1.1.
In this section we describe how to define complex numbers in VBA and how to construct functions for operations on complex numbers. Note that it is possible to use the built-in complex number functions in Excel directly, without having to construct them in VBA. However, we will see in later chapters that using the built-in functions increases substantially the computation time required for convergence of option prices. Constructing complex numbers in VBA, therefore, makes computation of option prices more efficient. Moreover, it is sometimes preferable to have control over how certain operations on complex numbers are defined. There are other definitions of the square root of a complex number, for example, than that given by applying DeMoivre’s Theorem. Finally, learning how to construct complex numbers in VBA is a good learning exercise.
The Excel file Chapter1Complex contains VBA functions to define complex numbers and to perform operations on complex numbers. Each function returns the real part and the imaginary part of the resulting complex number. The first step is to construct a complex number in terms of its two parts. The function Set_cNum() defines a complex number with real and imaginary parts given by set_cNum.rp and set_cNum.ip, respectively.
Function Set_cNum(rPart, iPart) As cNum
Set_cNum.rP = rPart
Set_cNum.iP = iPart
End Function
The function cNumProd() multiplies two complex numbers cNum1 and cNum2, and returns the complex number cNumProd with real and imaginary parts cNumProd.rp and cNumProd.ip, respectively.
Function cNumProd(cNum1 As cNum, cNum2 As cNum) As cNum
cNumProd.rP = (cNum1.rP * cNum2.rP) - (cNum1.iP * cNum2.iP)
cNumProd.iP = (cNum1.rP * cNum2.iP) + (cNum1.iP * cNum2.rP)
End Function
Similarly, the functions cNumDiv(), cNumAdd(), and cNumSub() return the real and imaginary parts of a complex number obtained by, respectively, division, addition, and subtraction of two complex numbers, while the function cNumConj() returns the conjugate of a complex number.
The function cNumSqrt() returns the square root of a complex number:
Function cNumSqrt(cNum1 As cNum) As cNum
r = Sqr(cNum1.rP ^ 2 + cNum1.iP ^ 2)
y = Atn(cNum1.iP / cNum1.rP)
cNumSqrt.rP = Sqr(r) * Cos(y / 2)
cNumSqrt.iP = Sqr(r) * Sin(y / 2)
End Function
The functions cNumExp() and cNumLn() produce, respectively, the exponential of a complex number and the natural logarithm of a complex number using the VBA function Atn() for the inverse tan function (arctan).
Function cNumExp(cNum1 As cNum) As cNum
cNumExp.rP = Exp(cNum1.rP) * Cos(cNum1.iP)
cNumExp.iP = Exp(cNum1.rP) * Sin(cNum1.iP)
End Function
Function cNumLn(cNum1 As cNum) As cNum
r = (cNum1.rP^2 + cNum1.iP^2)^0.5
theta = Atn(cNum1.iP / cNum1.rP)
cNumLn.rP = Application.Ln(r)
cNumLn.iP = theta
End Function
Finally, the functions cNumReal() and cNumIm() return the real and imaginary parts of a complex number, respectively.
The Excel file Chapter1Complex illustrates how these functions work. The VBA function Complexop2() performs operations on two complex numbers:
Function Complexop2(rP1, iP1, rP2, iP2, operation)
Dim cNum1 As cNum, cNum2 As cNum, cNum3 As cNum
Dim output(2) As Double
cNum1 = setcnum(rP1, iP1)
cNum2 = setcnum(rP2, iP2)
Select Case operation
Case 1: cNum3 = cNumAdd(cNum1, cNum2) ' Addition
Case 2: cNum3 = cNumSub(cNum1, cNum2) ' Subtraction
Case 3: cNum3 = cNumProd(cNum1, cNum2) ' Multiplication
Case 4: cNum3 = cNumDiv(cNum1, cNum2) ' Division
End Select
output(1) = cNum3.rP
output(2) = cNum3.iP
complexop2 = output
End Function
The Complexop2() function requires five inputs, a real and imaginary part for each number, and the parameter corresponding to the operation being performed (1 through 4). Its output is an array of dimension two, containing the real and imaginary parts of the complex number. Figure 1.1 illustrates how this function works. To add the two numbers 11 + 3i and −3 + 4i, which appear in ranges C4:D4 and C5:D5 respectively, in cell C6 we type
FIGURE 1.1 Operations on Complex Numbers
and copy to cell D6, which produces the complex number 8 + 7i. Note that the output of the Complexop2() function is an array. The appendix to this book explains in detail how to output arrays from functions. Note also that the last argument of the function Complexop2() is cell F6, which contains the operation number (1) corresponding to addition.
Similarly, the function Complexop1() performs operations on a single complex number, in this example 4 + 5i. To obtain the complex conjugate, in cell C15 we type
and copy to cell D15 This is illustrated in the bottom part of Figure 1.1.
Complex numbers are abstract entities, but they are extremely useful because they can be used in algebraic calculations to produce solutions that are tangible. In particular, the option pricing models covered in this book require a probability density function for the logarithm of the stock price, X = log(S). From a theoretical standpoint, however, it is often easier to obtain the characteristic function φX(t) for log(S), given by
where
The probability density function for the logarithm of the stock price can then be obtained by inversion of φX(t):
One corollary of Levy’s inversion formula—an alternate inversion formula—is that the cumulative density function FX(x) = Pr(X < x) for the logarithm of the stock price can be obtained. The following expression is often used for the risk-neutral probability that a call option lies in-the-money:
where k = log(K) is the logarithm of the strike price K. Again, this formula requires evaluating an integral that contains 
.
In this section we present two algorithms for finding roots of functions, the Newton-Raphson method, and the bisection method. These will become important in later chapters that deal with Black-Scholes implied volatility. Since the Black-Scholes formula cannot be inverted to yield the volatility, finding implied volatility must be done numerically. For a given market price on an option, implied volatility is that volatility which, when plugged into the Black-Scholes formula, produces the same price as the market. Equivalently, implied volatility is that which produces a zero difference between the market price and the Black-Scholes price. Hence, finding implied volatility is essentially a root-finding problem.
The chief advantage of programming root-finding algorithms in VBA, rather than using the Goal Seek and Solver features included in Excel, is that a particular algorithm can be programmed for the problem at hand. For example, we will see in later chapters that the bisection algorithm is particularly well suited for finding implied volatility. There are at least four considerations that must be kept in mind when implementing root-finding algorithms. First, adequate starting values must be carefully chosen. This is particularly important in regions of highly functional variability and when there are multiple roots and local minima. If the function is highly variable, a starting value that is not close enough to the root might stray the algorithm away from a root. If there are multiple roots, the algorithm may yield only one root and not identify the others. If there are local minima, the algorithm may get stuck in a local minimum. In that case, it would yield the minimum as the best approximation to the root, without realizing that the true root lies outside the region of the minimum. Second, the tolerance must be specified. The tolerance is the difference between successive approximations to the root. In regions where the function is flat, a high number for tolerance can be used. In regions where the function is very steep, however, a very small number must be used for tolerance. This is because even small deviations from the true root can produce values for the function that are substantially different from zero. Third, the maximum number of iterations needs to be defined. If the number of iterations is too low, the algorithm may stop before the tolerance level is satisfied. If the number of iterations is too high and the algorithm is not converging to a root because of an inaccurate starting value, the algorithm may continue needlessly and waste computing time.
To summarize, while the built-in modules such as the Excel Solver or Goal Seek allows the user to specify starting values, tolerance, maximum number of iterations, and constraints, writing VBA functions to perform root finding sometimes allows flexibility that built-in modules do not. Furthermore, programming multivariate optimization algorithms in VBA, such as the Nelder-Mead covered later in this chapter, is easier if one is already familiar with programming single-variable algorithms. The root-finding methods outlined in this section can be found in Burden and Faires (2001) or Press et al. (2002).
This method is one of the oldest and most popular methods for finding roots of functions. It is based on a first-order Taylor series approximation about the root. To find a root x of a function f(x), defined as that x which produces f(x) = 0, select a starting value x0 as the initial guess to the root, and update the guess using the formula
for i = 0, 1, 2, …, and where f′(xi) denotes the first derivative of f(x) evaluated at xi. There are two methods to specify a stopping condition for this algorithm, when the difference between two successive approximations is less than the tolerance level ε, or when the slope of the function is sufficiently close to zero. The VBA code in this chapter uses the second condition, but the code can easily be adapted for the first condition.
The Excel file Chapter1Roots contains the VBA functions for implementing the root-finding algorithms presented in this section. The file contains two functions for implementing the Newton-Raphson method. The first function assumes that an analytic form for the derivative f′(xi) exists, while the second uses an approximation to the derivative. Both are illustrated with the simple function f(x) = x2 − 7x + 10, which has the derivative f′(x) = 2x − 7. These are defined as the VBA functions Fun1() and dFun1(), respectively.
Function Fun1(x)
Fun1 = x^2 - 7*x + 10
End Function
Function dFun1(x)
dFun1 = 2*x - 7
End Function
The function NewtRap() assumes that the derivative has an analytic form, so it uses the function Fun1() and its derivative dFun1() to find the root of Fun1. It requires as inputs the function, its derivative, and a starting value x_guess. The maximum number of iterations is set at 500, and the tolerance is set at 0.00001.
Function NewtRap(fname As String, dfname As String, x_guess)
Maxiter = 500
Eps = 0.00001
cur_x = x_guess
For i = 1 To Maxiter
fx = Run(fname, cur_x)
dx = Run(dfname, cur_x)
If (Abs(dx) < Eps) Then Exit For
cur_x = cur_x - (fx / dx)
Next i
NewtRap = cur_x
End Function
The function NewRapNum() does not require the derivative to be specified, only the function Fun1() and a starting value. At each step, it calculates an approximation to the derivative.
Function NewtRapNum(fname As String, x_guess)
Maxiter = 500
Eps = 0.000001
delta_x = 0.000000001
cur_x = x_guess
For i = 1 To Maxiter
fx = Run(fname, cur_x)
fx_delta_x = Run(fname, cur_x - delta_x)
dx = (fx - fx_delta_x) / delta_x
If (Abs(dx) < Eps) Then Exit For
cur_x = cur_x - (fx / dx)
Next i
NewtRapNum = cur_x
End Function
The function NewtRapNum() approximates the derivative at any point x by using the line segment joining the function at x and at x + dx, where dx is a small number set at 1 × 10−9. This is the familiar “rise over run” approximation to the slope, based on a first-order Taylor series expansion for f(x + dx) about x:
This approximation appears as the statement
dx = (fx - fx_delta_x) / delta_x
in the function NewtRapNum().
This method is well suited to problems for which the function is continuous on an interval [a, b] and for which the function is known to take a positive value on one endpoint and a negative value on the other endpoint. By the Intermediate Value Theorem, the interval will necessarily contain a root. A first guess for the root is the midpoint of the interval. The bisection algorithm proceeds by repeatedly dividing the subintervals of [a, b] in two, and at each step locating the half that contains the root. The function BisMet() requires as inputs the function for which a root must be found, and the endpoints a and b. The endpoints must be chosen so that the function assumes opposite signs at each, otherwise the algorithm may not converge.
Function BisMet(fname As String, a, b)
Eps = 0.000001
If (Run(fname, b) < Run(fname, a)) Then
tmp = b: b = a: a = tmp
End If
Do While (Run(fname, b) - Run(fname, a) > Eps)
midPt = (b + a) / 2
If Run(fname, midPt) < 0 Then
a = midPt
Else
b = midPt
End If
Loop
BisMet = (b + a) / 2
End Function
We will see in Chapters 4 and 10 that the bisection method is particularly well suited for finding implied volatilities extracted from option prices.
Figure 1.2 illustrates the Newton-Raphson method with an explicit derivative, the Newton-Raphson method with an approximation to the derivative, and the Bisection method. This spreadsheet appears in the Excel file Chapter1Roots. As before, we use the function f(x) = x2 − 7x + 10, coded by the VBA function Fun1(), with derivative f′(x) = 2x − 7, coded by the VBA function dFun1(). It is easy to see by inspection that this function has two roots, at x = 2 and at x = 5. We illustrate the methods with the first root.
FIGURE 1.2 Root-Finding Algorithms
The bisection method requires an interval with endpoints chosen so that the function takes on values opposite in sign at each endpoint. Hence, we choose the interval [1, 3] for the first root, which appears in cells E7:E8. In cell G7 we type
which yields the value 4 for the function evaluated at the point x = 1. Similarly, in cell G8 we obtain the value −2 for the function evaluated at x = 3.
Recall that the VBA function BisMet() requires three inputs, a function name enclosed in quotes, and the endpoints of the interval along the x-axis over which the function changes sign. To invoke the bisection method, therefore, in cell C7 we type
which produces the root x = 2.
To invoke the two Newton-Raphson methods, we choose x0 = 1 as the starting value for the root x = 2. The VBA function NewtRap() uses an explicit form for the derivative and requires three inputs, the function name and the derivative name, each enclosed in quotes, and the starting value. Hence, in cell C8 we type
and obtain the root x = 2.
The VBA function NewRapNum() does not use an explicit form for the derivative, so it requires as inputs only the function name and a starting value. Hence, in cell C9 we type
and again obtain the root x = 2. The other root x = 5 is obtained similarly, using the interval [4,7] for the bisection algorithm and the starting value x0 = 4.
This example illustrates that proper selection of starting values and intervals is crucial, especially when multiple roots are involved. With the bisection method, an interval over which the function changes sign must be found for every root. Sometimes no such interval can be found, as is the case for the function f(x) = x2, which has a root at x = 0 but which never takes on negative values. In that case, the bisection method cannot be used.
With Newton’s method, it is important to select starting values close enough to every root that must be found. If not, the method might focus on one root only, and multiple roots may never be identified. Unfortunately, there is no method to properly identify appropriate starting values, so these are usually found by trail and error. In the case of a single variable, covered in this chapter, this is relatively straightforward and can be accomplished by dividing the x-axis into a series of starting values, and invoking Newton’s method at every starting value. In the multidimensional case, however, more complicated grid-search algorithms must be used.
In this section we present VBA code to perform multiple regression analysis under ordinary least squares (OLS) and weighted least squares (WLS). While functions to perform multiple regression are built into Excel, it is sometimes preferable to write VBA code to run regression, rather than relying on the built-in functions. First, Excel estimates parameters by OLS only, according to which each observation receives equal weight. If the analyst feels more weight should be given to certain observations, and less to others, then estimating parameters by WLS is preferable to OLS. Second, it is straightforward to obtain the entire covariance matrix of parameter estimates with VBA, rather than just its diagonal elements, which are the variance of each parameter estimate. Third, it is easy to obtain the “hat” matrix, whose diagonal elements can help identify the relative influence of each observation on parameter estimates. Finally, changing the contents of one cell automatically updates the results when a function is used to implement regression. This is not the case when OLS is implemented with the built-in regression routine in Excel.
To summarize, writing VBA code to perform regression is more flexible and allows the analyst to have access to more diagnostic tools than relying on the regression functions built into Excel. The OLS and WLS methods are explained in textbooks such as those by Davidson and MacKinnon (1993) and Neter et al. (1996).
Suppose we specify that the dependent variable Y is related to k − 1 independent variables X1, X2, …, Xk−1 and an intercept β0 in the linear form
where
The ordinary least-squares estimate of the parameters is given by the well-known formula
where
Once the OLS parameter estimates are obtained, we can obtain the fitted values as the vector 
of dimension n. If a model with no intercept is desired, then
and the design matrix excludes the vector ι containing ones, resulting in a design matrix of dimension n × (k − 1), and the OLS parameters are estimated by (1.4) as before.
It is very convenient to break down the total sum of squares (SSTO), defined as the variability of the dependent variable about its mean, in terms of the error sum of squares (SSE) and the regression sum of squares (SSR). Using algebra, it is straightforward to show that
With these quantities we can obtain several common definitions. An estimate of the variance σ2 is given by the mean square error (MSE)
while an estimate of the standard deviation is given by 
. The coefficient of multiple determination, R2, is given by the proportion of SSTO explained by the model
The R2 coefficient is the proportion of variability in the dependent variable that can be attributed to the linear model. The rest of the variability cannot be attributed to the model and is therefore pure unexplained error. One shortcoming of R2 is that it always increases, and never decreases, when additional independent variables are included in the model, regardless of whether or not the added variables have any explanatory power. The adjusted R2 incorporates a penalty for additional variables, and is given by
An estimate of the (k × k) covariance matrix of the parameter estimates is given by
The t-statistics for each regression coefficient are given by
where
Each t-statistic is distributed as a t random variable with n − k degrees of freedom, and can be used to perform a two-tailed test that each regression coefficient is zero. The two-tailed p-value is used to assess statistical significance of the regression coefficient. A small p-value corresponds to a coefficient that is significantly different from zero, whereas a large p-value denotes a coefficient that is statistically indistinguishable from zero. Usually 
is taken as the cut-off value to determine significance of each coefficient, corresponding to a significance level of 5 percent.
When an intercept term is included in the model, it can be shown by algebra that the condition 
always holds. When no intercept term is included, however, this condition may or may not hold. It is possible to obtain negative values of 
, especially if SSR is very small relative to SSTO.
Ordinary least squares attribute equal weight to each observation. In certain instances, however, the analyst may wish to assign more weight to some observations, and less weight to others. In this case, WLS is preferable to OLS. Selecting the weights w1, w2, …, wn, however, is arbitrary. One popular choice is to choose the weights as the inverse of each observation. Observations with a large variance receive little weight, while observations with a small variance get large weight.
Obtaining parameter estimates and associated statistics of (1.3) under WLS is straightforward. Define W as a diagonal matrix of dimension n containing the weights, so that W = diag[w1, …, wn]. Parameter estimates under WLS are given by
while SSTO, SSE, and SSR are given by
The (k × k) covariance matrix is given by
where MSE is given by (1.5), but using the definition of SSE given in (1.11), and the t-statistics are given by (1.9), but using the standard errors obtained as the square root of the diagonal elements of (1.12). The coefficients R2 and 
are given by (1.6) and (1.7) respectively, but using the sums of squares given in (1.11).
Under WLS, R2 does not have the convenient interpretation that it does under OLS. Hence, R2 and 
must be used with caution when these are obtained by WLS. Finally, we note that WLS is a special case of generalized least squares (GLS), according to which the matrix W is not a diagonal matrix but rather a matrix of general form. Under GLS, for example, the independence of the error terms can be relaxed to allow for different dependence structures between the errors.
In this section we present the VBA code for implementing OLS and WLS. We illustrate this with a simple example involving two explanatory variables, a vector of weights, and an intercept. The Excel file Chapter1WLS contains VBA code to implement WLS, and OLS as a special case. The function Diag() creates a diagonal matrix using a vector of weights as inputs:
Function Diag(W) As Variant
Dim n, i, j, k As Integer
Dim temp As Variant
n = W.Count
ReDim temp(n, n)
For i = 1 To n
For j = 1 To n
If j = i Then temp(i, j) = W(i) Else temp(i, j) = 0
Next j
Next i
Diag = temp
End Function
The function WLSregress() performs weighted least squares, and requires as inputs a vector y of observations for the dependent variable; a matrix X for the independent variables (which will contain ones in the first column if an intercept is desired); and a vector W of weights. This function produces WLS estimates of the regression parameters, given by (1.10). It is useful when only parameter estimates are required.
Function WLSregress(y As Variant, X As Variant, W As Variant) As Variant
Wmat = Diag(W)
n = W.Count
Dim Xtrans, Xw, XwX, XwXinv, Xwy As Variant
Dim m1, m2, m3, m4 As Variant
Dim output() As Variant
Xtrans = Application.Transpose(X)
Xw = Application.MMult(Xtrans, Wmat)
XwX = Application.MMult(Xw, X)
XwXinv = Application.MInverse(XwX)
Xwy = Application.MMult(Xw, y)
b = Application.MMult(XwXinv, Xwy)
k = Application.Count(b)
ReDim output(k) As Variant
For bcnt = 1 To k
output(bcnt) = b(bcnt, 1)
Next bcnt
WLSregress = Application.Transpose(output)
End Function
Note that the first part of the function creates a diagonal matrix Wmat of dimension n using the function Diag(), while the second part computes the WLS parameter estimates given by (1.10).
The second VBA function, WLSstats(), provides a more thorough WLS analysis and is useful when both parameter estimates and associated statistics are needed. It computes WLS parameter estimates, the standard error and t-statistic of each parameter estimate, its corresponding p-value, the R2 and 
coefficients, and 
, the estimate of the error standard deviation σ.
Function WLSstats(y As Variant, X As Variant, W As Variant) As Variant
Wmat = diag(W)
n = W.Count
Dim Xtrans, Xw, XwX, XwXinv, Xwy As Variant
Dim btemp As Variant
Dim output() As Variant, r(), se(), t(), pval() As Double
Xtrans = Application.Transpose(X)
Xw = Application.MMult(Xtrans, Wmat)
XwX = Application.MMult(Xw, X)
XwXinv = Application.MInverse(XwX)
Xwy = Application.MMult(Xw, y)
b = Application.MMult(XwXinv, Xwy)
n = Application.Count(y)
k = Application.Count(b)
ReDim output(k, 7) As Variant, r2(n), ss(n), t(k), se(k), pval(k)
As Double
yhat = Application.MMult(X, b)
For ncnt = 1 To n
r2(ncnt) = Wmat(ncnt, ncnt) * (y(ncnt) - yhat(ncnt, 1)) ^ 2
ss(ncnt) = Wmat(ncnt, ncnt) * (y(ncnt) - Application.Average(y)) ^ 2
Next ncnt
sse = Application.Sum(r2): mse = sse / (n - k)
rmse = Sqr(mse): sst = Application.Sum(ss)
rsquared = 1 - sse / sst
adj_rsquared = 1 - (n - 1) / (n - k) * (1 - rsquared)
For kcnt = 1 To k
se(kcnt) = (XwXinv(kcnt, kcnt) * mse) ^ 0.5
t(kcnt) = b(kcnt, 1) / se(kcnt)
pval(kcnt) = Application.TDist(Abs(t(kcnt)), n - k, 2)
output(kcnt, 1) = b(kcnt, 1)
output(kcnt, 2) = se(kcnt)
output(kcnt, 3) = t(kcnt)
output(kcnt, 4) = pval(kcnt)
Next kcnt
output(1, 5) = rsquared
output(1, 6) = adj_rsquared
output(1, 7) = rmse
For i = 2 To k
For j = 5 To 7
output(i, j) = “ ”
Next j
Next i
WLSstats = output
End Function
As in the previous VBA function WLSregress(), the function WLSstats() first creates a diagonal matrix Wmat using the vector of weights specified in the input argument W. It then creates estimated regression coefficients under WLS and the associated statistics. The last loop ensures that cells with no output remain empty.
The WLSregress() and WLSstats() functions are illustrated in Figure 1.3, using n = 16 observations. Both functions require as inputs the vectors containing values of the dependent variable, of the independent variables, and of the weights. The first eight observations have been assigned a weight of 1, while the remaining eight observations have been assigned a weight of 2.
FIGURE 1.3 Weighted Least Squares
The first VBA function, WLSregress(), produces estimated weighted coefficients only. The dependent variable is contained in cells C4:C19, the independent variables (including the intercept) in cells D4:F19, and the weights in cells B4:B19. Hence, in cell G9 we type
and copy down to cells G10:G11, which produces the three regression coefficients estimated by WLS, namely 
.
The second VBA function, WLSstats(), produces more detailed output. It requires the same inputs as the previous function, so in cell G4 we type
and copy to the range G4:M6, which produces the three estimated regression coefficients, their standard errors, t-statistics and p-values, the R2 and 
coefficients and 
.
It is easy to generalize the two WLS functions for more than two independent variables. The Excel file Chapter1WLSbig contains an example of the WLSregress() and WLSstats() functions using five independent variables and an intercept. This is illustrated in Figure 1.4.
FIGURE 1.4 Weighted Least Squares with Five Independent Variables
The range of independent variables (including intercept) is now contained in cells D4:I19. Hence, to obtain the WLS regression coefficients in cell B23 we type
and copy to cells B23:B28. To obtain detailed statistics, in cell D23 we type
and copy to cells D23:J28.
Suppose that OLS estimates of the regression coefficients are needed, instead of WLS estimates. Then cells B4:B19 are filled with ones, corresponding to equal weighting for all observations and a weighing matrix that is the identity matrix, so that the WLS estimates (1.10) will reduce to the OLS estimates (1.4). This is illustrated in Figure 1.5 with the Excel file Chapter1WLSbig.
FIGURE 1.5 Ordinary Least Squares with Five Independent Variables
Using the built-in regression data analysis module in Excel, it is easy to verify that Figure 1.5 produces the correct regression coefficients and associated statistics.
Finally, suppose that the intercept β0 is to be excluded from the model (1.3). In that case, the column corresponding to the intercept is excluded from the set of independent variables. This is illustrated in Figure 1.6.
FIGURE 1.6 Ordinary Least Squares with No Intercept
The range of independent variables is contained in cells E4:I19. Hence, in cell D23 we type
and copy to the range D23:J27, which produces the regression coefficients and associated statistics.
The methods to find roots of functions described earlier in this chapter are applicable when the root or minimum value of a single-variable function needs to be found. In many option pricing formulas, however, the minimum or maximum value of a function of two or more variables must be found. The Nelder-Mead algorithm is a powerful and popular method to find roots of multivariate functions. It is easy to implement, and it converges very quickly regardless of which starting values are used. Many mathematical and engineering packages, such as Matlab™ for example, use the Nelder-Mead algorithm in their optimization routines. We follow the description of the algorithm presented in Lagarias et al. (1999) for finding the minimum value of a multivariate function. Finding the maximum of a function can be done by changing the sign of the function and finding the minimum of the changed function.
For a function f(x) of n variables, the algorithm requires n + 1 starting values in x. Arrange these n + 1 starting values in increasing value for f(x), so that x1, x2, …, xn+1 are such that
where fk ≡ f(xk) and xi ∈ ℜn (i = 1, 2, …, n + 1). The best of these vectors is x1 since it produces the smallest value of f(x), and the worst is xn+1 since it produces the largest value. The remaining vectors lie in the middle. At each iteration step, the best values x1, …, xn are retained and the worst xn+1 is replaced according to the following rules:
where 
is the mean of the best n points and evaluate fr = f(xr). If 
then xn+1 is replaced with xr, the n + 1 points x1, …, xn, xr are reordered according to the value of the function as in (1.13), which produces another set of ordered points x1, x2, …, xn+1. The next iteration is initiated on the new worst point xn+1. Otherwise, proceed to the next rule.
and the value of the function fe = f(xe). If fe < fr then replace xn+1 with xe, reorder the points and initiate the next iteration. Otherwise, proceed to the next rule.
compute the outside contraction point 
and the value foc = f(xoc). If 
then replace xn+1 with xoc, reorder the points, and initiate the next iteration. Otherwise, proceed to rule 5 and perform a shrink step.
compute the inside contraction point 
and the value fic = f(xic). If fic < fn+1 replace xn+1 with xic, reorder the points and initiate the next iteration. Otherwise, proceed to rule 5.
for i = 2, …, n + 1. The new unordered points are the n + 1 points x1, v2, v3, …, vn+1. Reorder these points and initiate the next iteration. The Excel file Chapter1NM contains the VBA function NelderMead() for implementing this algorithm. The function uses a bubble sort algorithm to sort the values of the function in accordance with (1.13). This algorithm is implemented with the function BubSortRows().
Function BubSortRows(passVec)
Dim tmpVec() As Double, temp() As Double
uVec = passVec
rownum = UBound(uVec, 1)
colnum = UBound(uVec, 2)
ReDim tmpVec(rownum, colnum) As Double
ReDim temp(colnum) As Double
For i = rownum - 1 To 1 Step -1
For j = 1 To i
If (uVec(j, 1) > uVec(j + 1, 1)) Then
For k = 1 To colnum
temp(k) = uVec(j + 1, k)
uVec(j + 1, k) = uVec(j, k)
uVec(j, k) = temp(k)
Next k
End If
Next j
Next i
BubSortRows = uVec
End Function
In this chapter the Nelder-Mead algorithm is illustrated using the bivariate function f(x1, x2) = f(x, y) defined as
and the function of three variables g(x1, x2, x3) = g(x, y, z) defined as
These are coded as the VBA functions Fun1() and Fun2(), respectively.
Function Fun1(params)
x = params(1)
y = params(2)
Fun1 = x ^ 2 - 4 * x + y ^ 2 - y - x * y
End Function
Function Fun2(params)
x = params(1)
y = params(2)
z = params(3)
Fun2 = (x - 10) ^ 2 + (y + 10) ^ 2 + (z - 2) ^ 2
End Function
The function NelderMead() requires as inputs only the name of the VBA function for which a minimum is to be found, and a set of starting values.
Function NelderMead(fname As String, startParams)
Dim resMatrix() As Double
Dim x1() As Double, xn() As Double, xw() As Double, xbar()
As Double, xr() As Double, xe() As Double, xc() As Double,
xcc() As Double
Dim funRes() As Double, passParams() As Double
MAXFUN = 1000
TOL = 0.0000000001
rho = 1
Xi = 2
gam = 0.5
sigma = 0.5
paramnum = Application.Count(startParams)
ReDim resmat(paramnum + 1, paramnum + 1) As Double
ReDim x1(paramnum) As Double, xn(paramnum) As Double,
xw(paramnum) As Double, xbar(paramnum) As Double,
xr(paramnum) As Double, xe(paramnum) As Double,
xc(paramnum) As Double, xcc(paramnum) As Double
ReDim funRes(paramnum + 1) As Double, passParams(paramnum)
For i = 1 To paramnum
resmat(1, i + 1) = startParams(i)
Next i
resmat(1, 1) = Run(fname, startParams)
For j = 1 To paramnum
For i = 1 To paramnum
If (i = j) Then
If (startParams(i) = 0) Then
resmat(j + 1, i + 1) = 0.05
Else
resmat(j + 1, i + 1) = startParams(i) * 1.05
End If
Else
resmat(j + 1, i + 1) = startParams(i)
End If
passParams(i) = resmat(j + 1, i + 1)
Next i
resmat(j + 1, 1) = Run(fname, passParams)
Next j
For j = 1 To paramnum
For i = 1 To paramnum
If (i = j) Then
resmat(j + 1, i + 1) = startParams(i) * 1.05
Else
resmat(j + 1, i + 1) = startParams(i)
End If
passParams(i) = resmat(j + 1, i + 1)
Next i
resmat(j + 1, 1) = Run(fname, passParams)
Next j
For lnum = 1 To MAXFUN
resmat = BubSortRows(resmat)
If (Abs(resmat(1, 1) - resmat(paramnum + 1, 1)) < TOL) Then
Exit For
End If
f1 = resmat(1, 1)
For i = 1 To paramnum
x1(i) = resmat(1, i + 1)
Next i
fn = resmat(paramnum, 1)
For i = 1 To paramnum
xn(i) = resmat(paramnum, i + 1)
Next i
fw = resmat(paramnum + 1, 1)
For i = 1 To paramnum
xw(i) = resmat(paramnum + 1, i + 1)
Next i
For i = 1 To paramnum
xbar(i) = 0
For j = 1 To paramnum
xbar(i) = xbar(i) + resmat(j, i + 1)
Next j
xbar(i) = xbar(i) / paramnum
Next i
For i = 1 To paramnum
xr(i) = xbar(i) + rho * (xbar(i) - xw(i))
Next i
fr = Run(fname, xr)
shrink = 0
If ((fr >= f1) And (fr < fn)) Then
newpoint = xr
newf = fr
ElseIf (fr < f1) Then
For i = 1 To paramnum
xe(i) = xbar(i) + Xi * (xr(i) - xbar(i))
Next i
fe = Run(fname, xe)
If (fe < fr) Then
newpoint = xe
newf = fe
Else
newpoint = xr
newf = fr
End If
ElseIf (fr >= fn) Then
If ((fr >= fn) And (fr < fw)) Then
For i = 1 To paramnum
xc(i) = xbar(i) + gam * (xr(i) - xbar(i))
Next i
fc = Run(fname, xc)
If (fc <= fr) Then
newpoint = xc
newf = fc
Else
shrink = 1
End If
Else
For i = 1 To paramnum
xcc(i) = xbar(i) - gam * (xbar(i) - xw(i))
Next i
fcc = Run(fname, xcc)
If (fcc < fw) Then
newpoint = xcc
newf = fcc
Else
shrink = 1
End If
End If
End If
If (shrink = 1) Then
For scnt = 2 To paramnum + 1
For i = 1 To paramnum
resmat(scnt, i + 1) = x1(i) + sigma
* (resmat(scnt, i + 1) - x1(1))
passParams(i) = resmat(scnt, i + 1)
Next i
resmat(scnt, 1) = Run(fname, passParams)
Next scnt
Else
For i = 1 To paramnum
resmat(paramnum + 1, i + 1) = newpoint(i)
Next i
resmat(paramnum + 1, 1) = newf
End If
Next lnum
If (lnum = MAXFUN + 1) Then
MsgBox “Maximum Iteration (“ & MAXFUN & “) exceeded”
End If
resmat = BubSortRows(resmat)
For i = 1 To paramnum + 1
funRes(i) = resmat(1, i)
Next i
funRes(1) = funRes(1)
NelderMead = Application.Transpose(funRes)
End Function
In this function, the maximum number of iterations (MAXFUN) has been set to 1,000, and the tolerance (TOL) between the best value and worst values x1, and xn+1, respectively, to 10−10. Increasing MAXFUN and decreasing TOL will lead to more accurate results, but will also increase the computing time.
Figure 1.7 illustrates the NelderMead() function on the functions f and g defined in (1.14) and (1.15). It is easy to verify that f takes on its minimum value of −7 at (x, y) = (3, 2) and that g takes on its minimum value of 0 at the point (x, y, z) = (10, − 10, 2). We choose large starting values of 10,000 to illustrate the convergence of the Nelder-Mead algorithm.
FIGURE 1.7 Nelder-Mead Algorithm
The starting values for f(x, y) appear in cells F7:F8, and the final values in cells E7:E8. The value of f(3, 2) = −7 appears in cell E6. The starting values for g(x, y, z) appear in cells F14:F16, the final values in cells E14:E16 and the value of g(10, − 10, 2) = 0 in cell E13. Even with starting values that are substantially different than the optimized values, the Nelder-Mead algorithm is able to locate the minimum of each function.
All of the option pricing models presented in this book rely on parameters as input values to the models. Estimates of these parameters must therefore be found. One popular method to estimate parameter values is the method of maximum likelihood. This method is explained in textbooks on mathematical statistics and econometrics, such as that by Davidson and MacKinnon (1993). Maximum likelihood estimation requires a sample of n identically and independently distributed observations x1, x2, …, xn, assumed to originate from some parametric family. In that case, the joint probability density function of the sample (the likelihood function) is factored into a product of marginal functions, and the maximum likelihood estimators (MLEs) of the parameters are those values that maximize the likelihood.
In some cases the MLEs are available in closed form. For example, it is well known that if x1, x2, …, xn are from the normal distribution with parameters μ for the mean and σ2 for the variance, then Davidson and MacKinnon (1993) explain that the Type II maximum likelihood estimators of μ and σ2 can be found by differentiating the likelihood function. These have the closed form 
and 
respectively, where both summations run from i = 1 to n.
For most of the option pricing models covered in this book, however, MLEs of parameters are not available in closed form. Hence, we need to apply an optimization routine such as the Nelder-Mead algorithm to find those values of the parameters that maximize the likelihood. If x1, x2, …, xn are independently and identically distributed according to some distribution f(xi;θ), where θ denotes a vector of parameters, then the Type I MLE of θ is given by
In Chapter 6 we apply the Nelder-Mead algorithm to expressions such as (1.16), to find MLEs of parameters for generalized autoregressive conditional heteroskedastic (GARCH) models of volatility.
Sometimes functional values ai = f(xi) are available only on a discrete set of points xi (i = 1, …, n), even though values at all points are needed. In Chapter 10, for example, we will see that implied volatilities are usually calculated at strike prices obtained at intervals of $5, but that it is useful to have implied volatilities at all strike prices. Interpolation refers to a wide class of methods that allow functional values to be joined together in a piecewise fashion, so that a continuous graph can be obtained. Starting with a set of points 
, the cubic spline interpolating function is a series of cubic polynomials defined as
Each Sj(x) is defined for values of x in the subinterval [xj, xj+1], where there are n − 1 such subintervals. Finding the coefficients bj, cj and dj is done by solving a linear system of equations which can be represented as Ax = b, where A is a tridiagonal matrix. Natural cubic splines are produced when the boundary conditions are set as S″(x1) = S″(xn) = 0, while clamped cubic splines, when S″(x1) = f′(x1) and S″(xn) = f′(xn). To obtain natural cubic splines, define hj = xj+1 − xj for j = 1, …, n − 1, and for j = 2, …, n − 1 set 
. Then initialize ℓ1 = ℓn = 1, μ1 = μn = 0, z1 = zn = 0, and for j = 2, …, n − 1 define ℓj = 2(xj+1 − xj−1) − hj−1μj−1, μj = hj/ℓj, and zj = (αj − hj−1zj−1)/ℓj. Finally, working backward, the coefficients are obtained as
for j = n − 1, n − 2, …, 1. The Excel file Chapter1NatSpline contains the VBA function NSpline() for implementing natural cubic splines. The NSpline() function requires as inputs the points 
and the functional values 
. It returns S(x) as defined in (1.17) for any value of x.
Function NSpline(s, x, a)
n = Application.Count(x)
Dim h() As Double, alpha() As Double
ReDim h(n - 1) As Double, alpha(n - 1) As Double
For i = 1 To n - 1
h(i) = x(i + 1) - x(i)
Next i
For i = 2 To n - 1
alpha(i) = 3 / h(i) * (a(i + 1) - a(i)) - 3 / h(i - 1)
* (a(i) - a(i - 1))
Next i
Dim l() As Double, u() As Double, z() As Double, c() As Double, b() As Double, d() As Double
ReDim l(n) As Double, u(n) As Double, z(n) As Double, c(n) As Double, b(n) As Double, d(n) As Double
l(1) = 1: u(1) = 0: z(1) = 0
l(n) = 1: z(n) = 0: c(n) = 0
For i = 2 To n - 1
l(i) = 2 * (x(i + 1) - x(i - 1)) - h(i - 1) * u(i - 1)
u(i) = h(i) / l(i)
z(i) = (alpha(i) - h(i - 1) * z(i - 1)) / l(i)
Next i
For i = n - 1 To 1 Step -1
c(i) = z(i) - u(i) * c(i + 1)
b(i) = (a(i + 1) - a(i)) / h(i) - h(i) * (c(i + 1) + 2 * c(i)) / 3
d(i) = (c(i + 1) - c(i)) / 3 / h(i)
Next i
For i = 1 To n - 1
If (x(i) <= s) And (s <= x(i + 1)) Then
NSpline = a(i) + b(i) * (s - x(i))
+ c(i) * (s - x(i)) ^ 2 + d(i) * (s - x(i)) ^ 3
End If
Next i
End Function
Figure 1.8 illustrates this function, with x-values xj contained in cells B5:B13, and functional values aj in cells C5:C13.
FIGURE 1.8 Natural Cubic Splines
To obtain the value of (1.17) when x = 2.6 in cell E20, for example, in cell F20 we type
which produces S(2.6) = 0.7357.
To implement clamped cubic splines, a similar algorithm is used, except that α1 = 3(a2 − a1)/h1 − 3f′1, αn = 3f′n − 3(an − an−1)/hn−1, where f′1 = f′(x1) and f′n = f′(xn) need to be specified, and the intermediate coefficients are initialized as ℓ1 = 2h1, μ1 = 0.5, z1 = α1/ℓ1, ℓn = hn−1(2 − μn−1), zn = (αn − hn−1zn−1)/ℓn, and cn = zn. The values ℓj, μj, and zj for j = 2, …, n − 1, as well as the coefficients cj, bj, and dj for j = n − 1, n − 2, …, 1 are obtained as before. The Excel file Chapter1ClampSpline contains the VBA function CSpline() for implementing the clamped spline algorithm. It requires as additional inputs values f′1 and f′n for the derivatives, via the variables fp1 and fpN.
Function CSpline(s, x, a, fp1, fpN)
n = Application.Count(x)
Dim h() As Double, alpha() As Double
ReDim h(n - 1) As Double, alpha(n) As Double
For i = 1 To n - 1
h(i) = x(i + 1) - x(i)
Next i
alpha(1) = 3 * (a(2) - a(1)) / h(1) - 3 * fp1
alpha(n) = 3 * fpN - 3 * (a(n) - a(n - 1)) / h(n - 1)
For i = 2 To n - 1
alpha(i) = 3 / h(i) * (a(i + 1) - a(i)) - 3 / h(i - 1)
* (a(i) - a(i - 1))
Next i
Dim l() As Double, u() As Double, z() As Double, c() As Double, b() As Double, d() As Double
ReDim l(n) As Double, u(n) As Double, z(n) As Double, c(n) As Double, b(n) As Double, d(n) As Double
l(1) = 2 * h(1): u(1) = 0.5: z(1) = alpha(1) / l(1)
For i = 2 To n - 1
l(i) = 2 * (x(i + 1) - x(i - 1)) - h(i - 1) * u(i - 1)
u(i) = h(i) / l(i)
z(i) = (alpha(i) - h(i - 1) * z(i - 1)) / l(i)
Next i
l(n) = h(n - 1) * (2 - u(n - 1))
z(n) = (alpha(n) - h(n - 1) * z(n - 1)) / l(n)
c(n) = z(n)
For i = n - 1 To 1 Step -1
c(i) = z(i) - u(i) * c(i + 1)
b(i) = (a(i + 1) - a(i)) / h(i) - h(i) * (c(i + 1) + 2 * c(i)) / 3
d(i) = (c(i + 1) - c(i)) / 3 / h(i)
Next i
For i = 1 To n - 1
If (x(i) <= s) And (s <= x(i + 1)) Then
CSpline = a(i) + b(i) * (s - x(i))
+ c(i) * (s - x(i)) ^ 2 + d(i) * (s - x(i)) ^ 3
End If
Next i
End Function
Figure 1.9 illustrates this function, using the same data in Figure 1.8, and with f′1 = f′n = 0.10 in cells I3 and I4.
FIGURE 1.9 Clamped Cubic Splines
To obtain the value of (1.17) when x = 2.6, in cell F20 we type
which produces S(2.6) = 0.7317. Note that clamped cubic splines are slightly more accurate, since they incorporate information about the function, in terms of the first derivatives of the function at the endpoints x1 and xn. In this example, however, the values of 0.10 are arbitrary so there is no advantage gained by using clamped cubic splines over natural cubic splines. For a detailed explanation of cubic splines, see Burden and Faires (2001), or Press et al. (2002).
In this chapter we illustrate how VBA can be used to perform operations on complex numbers, to implement root-finding algorithms and the minimization of multivariate functions, to perform estimation by ordinary least squares and weighted least squares, to find maximum likelihood estimators, and to implement cubic splines. While Excel provides built-in functions to perform operations on complex numbers, to implement root-finding methods, and to perform regression, there are advantages to programming these in VBA rather than relying on the built-in functions. In the case of root-finding methods, users are not restricted to the methods built into the Goal Seek or Solver Excel modules but can select from a wide variety of root-finding algorithms. Numerical analysis textbooks, such as that by Burden and Faires (2001), present many of the root finding algorithms currently available, and explain cubic splines in detail. Implementing regression with VBA also has its advantages. The most notable is that weighted least squares estimation can be implemented so that users are not restricted to using only ordinary least squares. As well, users can program many additional features, such as diagnostic statistics or the covariance matrix, that are not available in Excel regression module built into Excel. One additional advantage is that regression results are updated automatically every time the user changes the values of a data cell, without having to rerun the regression.
This section provides exercises to reinforce the mathematical concepts introduced in this chapter. Solutions are contained in the Excel file Chapter1Exercises.
1.1 Write VBA functions to obtain the cosine of a complex number z = x + iy, whose real and imaginary parts are given by cos (x) cosh (y) and − sin (x) sinh (y) respectively, and the sine of a complex number, whose real and imaginary parts are sin (x) cosh (y) and − cos (x) sinh (y) respectively. The hyperbolic sine and hyperbolic cosine of a real number x are defined as sinh (x) = (ex − e−x)/2 and cosh (x) = (ex + e−x)/2.
1.2 Use DeMoivre’s Theorem (1.1) to find an expression for the power of a complex number, and write a VBA function to implement it.
1.3 Find an expression for a complex number raised to the power of another complex number. That is, find an expression for wz, where w and z are both complex numbers. This requires the angle formulas
1.4 The hat matrix H is very useful for identifying outliers in data. Its diagonal entries, called leverage values, are used to gauge the relative influence of each data point and to construct studentized residuals and deleted residuals. The hat matrix of dimension (n × n) and is given by
where I = identity matrix of dimension n × n, and
X = design matrix of dimension (n × k).
Write VBA code to create the hat matrix. In your code, output the diagonal elements of the hat matrix, hii (i = 1, 2, …, n) to an array in the spreadsheet. Use the data for two independent variables contained in the spreadsheet Chapter1WLS.
1.5 The global F-test of a regression model under OLS is a test of the null hypothesis that all regression coefficients are zero simultaneously. The F-statistic is defined as
where SSR is the regression sum of squares, SSE is the error sum of squares, n is the number of observations, and k is the number of parameters (including the intercept). Under the null hypothesis, F* is distributed as an F random variable with k − 1 and n − k degrees of freedom. Write VBA code to append the Excel file Chapter1WLS to include the global F-test. Include in your output a value of the F-statistic, and the p-value of the test. Use the Excel function Fdist() to find this p-value.
1.1 Define two new VBA functions cNumCos() and cNumSin() and add these at the bottom of the VBA code in the Excel file Chapter1Complex.
Function cNumCos(cNum1 As cNum) As cNum
cNumCos.rP = Cos(cNum1.rP) * Application.Cosh(cNum1.iP)
cNumCos.iP = -Sin(cNum1.rP) * Application.Sinh(cNum1.iP)
End Function
Function cNumSin(cNum1 As cNum) As cNum
cNumSin.rP = Sin(cNum1.rP) * Application.Cosh(cNum1.iP)
cNumSin.iP = -Cos(cNum1.rP) * Application.Sinh(cNum1.iP)
End Function
Modify the VBA function Complexop1() to include two additional cases for the Select Case statement
Case 5: cNum3 = cNumCos(cNum) ' Cosine
Case 6: cNum3 = cNumSin(cNum) ' Sine
The cosine and sine of a complex number can be obtained by invoking the VBA function Complexop1() described in Section 1.1.2. This is illustrated in Figure 1.10.
FIGURE 1.10 Solution to Exercise 1.1
1.2 Recall that a complex number w can be written w = r cos (y) + ir sin (y), where y = arctan(b/a) and 
. Using DeMoivre’s Theorem, w raised to the power n can be written as
Hence Re[wn] = rn cos (ny) and Im[wn] = rn sin (ny). The functions complexop3() and cNumPower() to obtain wn are presented below, and their use is illustrated in Figure 1.11.
FIGURE 1.11 Solution to Exercise 1.2
Function complexop3(rP, iP, n)
Dim cNum1 As cNum, cNum3 As cNum
Dim output(2) As Double
cNum1 = setcnum(rP, iP)
cNum3 = cNumPower(cNum1, n.Value)
output(1) = cNum3.rP
output(2) = cNum3.iP
complexop3 = output
End Function
Function cNumPower(cNum1 As cNum, n As Double) As cNum
r = Sqr(cNum1.rP ^ 2 + cNum1.iP ^ 2)
y = Atn(cNum1.iP / cNum1.rP)
cNumPower.rP = r ^ n * Cos(y * n)
cNumPower.iP = r ^ n * Sin(y * n)
End Function
For w = 4 + 5i and n = 2, the functions produce w2 = −9 + 40i. Note also that when n = 0.5 the functions produce the square root of a complex number.
1.3 We seek wz = (a + ib)c+id. Let 
, and θ = arctan(b/a), then
where rid = eid ln r. Now multiply and group trigonometric terms together
Finally, apply the angle formulas to arrive at
Hence, the real and imaginary parts of wz are
The VBA function cNumPowercNum() implements the exponentiation of a complex number by another complex number.
Function cNumPowercNum(cNum1 As cNum, cNum2 As cNum) As cNum
r = Sqr(cNum1.rP ^ 2 + cNum1.iP ^ 2)
y = Atn(cNum1.iP / cNum1.rP)
cNumPowercNum.rP = r ^ cNum2.rP * Exp(-cNum2.iP * y) *
Cos(cNum2.rP * y + cNum2.iP * Log(r))
cNumPowercNum.iP = r ^ cNum2.rP * Exp(-cNum2.iP * y) *
Sin(cNum2.rP * y + cNum2.iP * Log(r))
End Function
Function complexop4(rP1, iP1, rP2, iP2)
Dim cNum1 As cNum, cNum2 As cNum, cNum3 As cNum
Dim output(2) As Double
cNum1 = setcnum(rP1, iP1)
cNum2 = setcnum(rP2, iP2)
cNum3 = cNumPowercNum(cNum1, cNum2)
output(1) = cNum3.rP
output(2) = cNum3.iP
complexop4 = output
End Function
Use of this function is illustrated in Figure 1.12, using w = 4 + 5i and z = 2 + 3i. The results indicate that Re[wz] = 1.316 and Im[wz] = 2.458.
FIGURE 1.12 Solution to Exercise 1.3
1.4 The VBA function Hat() creates the hat matrix and outputs its diagonal entries to any cells in the spreadsheet:
Option Base 1
Function Hat(X As Variant) As Variant
Dim n, k As Integer
Dim output(), H() As Variant
n = X.Rows.Count
ReDim output(n, n), H(n, n) As Variant
Xt = Application.Transpose(X)
XtX_inv = Application.MInverse(Application.MMult(Xt, X))
XXtX_inv = Application.MMult(X, XtX_inv)
H = Application.MMult(XXtX_inv, Xt)
For k = 1 To n
output(k, 1) = H(k, k)
Next k
Hat = output
End Function
The results of this function appear in Figure 1.13. The function requires as input only the design matrix (including a column for the intercept), which appears in cells B7:D22. Hence, in cell F7 we type
FIGURE 1.13 Solution to Exercise 1.4
and copy the formula as an array to the remaining cells. It can be shown that the trace of H, the sum of its diagonal elements, is equal to k. Hence, in cell F24 we type “=sum(F7:F22)” and obtain k = 3. The diagonals of the hat matrix indicate that observation 9, with h9, 9 = 0.447412, has the largest influence, followed by observation 12 with h12, 12 = 0.287729, and so on.
1.5. To implement the global F-test, we must first increase the range for the output, which creates extra room for the F-statistic and its p-value. Hence, we change the first part of the ReDim statement from ReDim output(k,5) to ReDim output(k + 2, 5). Then we add the VBA code to create the F-statistic and its p-value using the Fdist() function in Excel, and two statements to output these quantities to the spreadsheet:
Dim F As Double
F = (sst - sse) / (k - 1) / mse
Fpval = Application.FDist(F, k - 1, n - k)
output(4, 5) = F
output(5, 5) = Fpval
For i = 4 To 5
For j = 1 To 4
output(i, j) = “ “
Next j
Next i
The last output statements ensure that zeros are not needlessly outputted to the spreadsheet. Figure 1.14 illustrates the global F-test using the data for two independent variables contained in the Excel file Chapter1WLS and presented earlier in this chapter (see Figure 1.3). These statements produce F* = 0.935329 and p = 0.417336. Using the built-in regression modules in Excel, or a statistical software package, it is easy to verify that these numbers are correct. The results of the F-test imply that we cannot reject the null hypothesis that the coefficients are all zero, so that the joint values of the coefficients are β0 = 0, β1 = 0, and β2 = 0.
FIGURE 1.14 Solution to Exercise 1.5
