Jan Hurt
Division of Financial & Insurance Mathematics
Faculty of Mathematics & Physics
CHARLES UNIVERSITY OF PRAGUE

Computer Algebra Systems (CAS) with Applications to Statistics and Finance

(Mathematica, Maple, Derive, ... )

4 pillars of Mathematica
1) Numerical calculations

2 + 2

N[π, 
      100]

2) Symbolic or algebraic calculations

Expand[(x + 2  y)^2]

3) Graphics

Plot3D[
        Sin[x + y^2], {x, -2, 2}, {y, -2, 2}];

4) Programming

fun[x_] := x^3Table[fun[i], {i, 1, 5}]

average[data_] := Apply[Plus, data]/Length[data]average[{1, 2, 3, 4, 5}]

\(TraditionalForm
    \`hilbert[n_] := Table[1/\((i + j - 1)\), {i, n}, {j, n}]\)\(TraditionalForm\`Inverse[hilbert[3]]
\)Short[Inverse[hilbert[30]][[1]]]

The reasons for using Mathematica in financial and insurance modelling 
(William Shaw, Quantitative Analysis Group, Nomura International plc, London):

(a) the fact that they have built-in routines for most basic common operations;
(b) the ability to do calculus with simple and complex functions, which can eliminate many of the headaches;
(c) the ability to compare analytical symbolic models directly with the results from a numerical algorithm;
(d) the ability to visualize the results, and literally see the errors if there are any.

Change-point (Linear-plateau), used in animal research

Needs["Statistics`Master`"]

linearplateau(x_, a_, b_, p_, c_) := 
      a - b p log((E^(((c - x)/p))) + 1)

({{({{linearplateau2(x_, a_, b_, p_, c_, k_) := 
                      a + k x - 
                        b p log((E^(((c - x)/p))) + 
                              1)}})}})

Plot[linearplateau[x, 10, 1,  .1, 10], {x, 0, 20}, 
        PlotRange -> All];

\(TraditionalForm\`\(Plot[
        linearplateau2[x, 10, 1,  .1, 13, 0.2], {x, 0, 20}, 
	
        PlotRange -> All];\)\)

distr1 = NormalDistribution[0, 1]; 
simul1 = 
            Table[Module[{}, part1 = Table[x + Random[distr1], {x, 10}]; 
                
part2 = Table[Random[distr1], {10}] + 10; 
                
dep = Join[part1, part2]; 
                
data1 = Transpose[{Range[20], dep}]; 
                
{data1, 
                  ListPlot[dep, PlotJoined -> True, 
                    PlotRange -> {-2, 12}, 
                    DisplayFunction -> Identity]}], {10}];

Show[
        Map[simul1[[#, 2]] &, Range[Length[simul1]]], 
        DisplayFunction -> $DisplayFunction];

results1= Map[ NonlinearFit[
Flatten[Map[Take[simul1[[#]],1]&,
Range[Length[simul1]]],1][[#]],
a-b*p*Log[1+Exp[(c-x)/p]],x,{a,b,p,c},
		 MaxIterations->30]&,
			Range[Length[simul1]]]

plots1 = 
        Map[Plot[#, {x, 0, 20}, PlotRange -> {-2, 12}, 
              DisplayFunction -> Identity] &, results1];

Show[
        Map[plots1[[#]] &, Range[Length[plots1]]], 
        DisplayFunction -> $DisplayFunction];

LAD (Least Absolute Deviation) Spaeth, H.: Algorithmen f. multivariable Ausgleichs-modelle p. 50

Beware of indices! Here m means the number of observations. Input: X(mxn) data matrix, y(mx1) vector of dependent variables. The procedure finds a(nx1) vector such that Xa-y is minimal in L1 norm.

Needs["LinearAlgebra`MatrixManipulation`"];
\(TraditionalForm\`L1Spaeth[X_, y_] := 
      Module[{m, n, im, nullsm, ones, x1, part1, 
						\ part2, c, 
          constr, opt, bint, aopt}, 

        m = \(Dimensions[X]\)[\([1]\)]; \  (*\ No . \ of\ observations\ *) 

        n = \(Dimensions[X]\)[\([2]\)]; \  (*\ No . \ of\ parameters\ *) 

        im = IdentityMatrix[m]; 
nullsm = ZeroMatrix[m]; 

        ones = Table[1, {n}]; 
x1 = X . ones; 

        part1 = Transpose[
            AppendColumns[\(-im\), im, Transpose[X], \(-{x1}\)]]; 

        part2 = Transpose[
            AppendColumns[im, \(-im\), \(-Transpose[X]\), {x1}]]; 

        matt = AppendColumns[part1, part2]; 

        c = Join[Table[1, {2*m}], Table[0, {n + 1}]]; 

        constr = Join[y, \(-y\)]; 

        optim = LinearProgramming[c, matt, constr]; 

        bint = Take[optim, \(-\((n + 1)\)\)]; 

        aopt = Take[bint, n] - Last[bint]\  // \ N
 (*\ 
          End\ of\ module*) ]\)

X = Transpose[{{1, 1, 1, 1, 1}, {3, 5, 1, 4, 2}}]; \(y1 = {2, 4, 1, 6, 4}; 
y2 = \(-{2, 4, 1, 6, 4}\); 
L1Spaeth[X, y1]\)L1Spaeth[X, y2]

\(InputForm
    \`p1 = ListPlot[
        Sort[{\(X\^T\)[[2]], y2}\^T], 
        PlotJoined -> True]; 

    p2 = Plot[\(- .75\)\ x -  .25, {x, 1, 5}]; 
Show[p1, p2]; \)

\(TraditionalForm
    \`p1 = ListPlot[Sort[Transpose[{\(Transpose[X]\)[\([2]\)], y2}]], 
        PlotJoined -> True, DisplayFunction -> Identity]; 

    p2 = Plot[\(-0.75\)*x - 0.25, {x, 1, 5}, 
        DisplayFunction -> Identity]; 

    Show[p1, p2, DisplayFunction -> $DisplayFunction]; \)

Example (Cauchy distribution, linear regression)

n = 50;  dist3 = CauchyDistribution[0, 1]; \(X = Transpose[{Table[1, {n}], Table[i, {i, n}]}]; 

      y = Table[0.5*i + Random[dist3], {i, n}]; 

      g1 = ListPlot[y, PlotJoined -> True, PlotRange -> All]; \)

s1 =   L1Spaeth[X, y]      s2 =   Fit[y, {1, x}, x]g2 = Plot[s1 . {1, x}, {x, 1, 50}, DisplayFunction -> Identity]
            ;
g3 = Plot[s2, {x, 1, 50}, DisplayFunction -> Identity];
Show[g2, g3, DisplayFunction -> $DisplayFunction];

Multidimensional Scaling

Scaling.nb 
Reference: Mardia, K. V. et al.: Multivariate Analysis. Academic Press 1979
(p. 397)
Version: January 8, 1999 
Copyright: Jan Hurt, 1999 

Notation: n number of points, p the dimension of the configuration space
makeconfiguration::usage="makeconfiguration[dit_,  pdim___] gives a configuration of points in a p-dimensional  Euclidean space where p is the number of positive eigenvalues of the respective matrix, or, if pdim given, p:=pdim. The construction uses the distances given in the form of lower diagonal list dit (without zeros on the diagonal), i.e., dit = {{(Subscript[d, 21])}, {(Subscript[d, 31]), (Subscript[d, 32])}, {(Subscript[d, 41]), (Subscript[d, 42]), (Subscript[d, 43])}, ... }. The configuration is in the form of (nxp) matrix with the centre of gravity (the average) 0."

newdistancematrix::usage="newdistancematrix[pot] returns the distance matrix of the points given in the (nxp) matrix pot. pot should be of the same form as the output of makeconfiguration."
		
newdistancediag::usage="newdistancediag[pot] returns the lower diagonal of the distance matrix of the points given in the (nxp) matrix pot. pot should be of the same form as the output of makeconfiguration."	

stress::usage="stress[distancediag, distancehatdiag] returns     
ψ = \(TraditionalForm\`∑\+\(i = 1\)\%n\(∑\+\(j = 1\)\%n\((d\_ij\)\)\) - \(TraditionalForm\`d\&^\)\(TraditionalForm\`\_ij\))^2, also called raw stress, where (Subscript[d, ij]) are the original distances and  \(TraditionalForm\`d\&^\)\(TraditionalForm\`\_ij\) are the "estimated" or reproduced distances from the configuration in p dimensions. Both of the parameters should be in the form of lower diagonal structures {{(Subscript[d, 21])}, {(Subscript[d, 31]), (Subscript[d, 32])}, {(Subscript[d, 41]), (Subscript[d, 42]), (Subscript[d, 43])}, ... } (without zeros) of the same dimension.

\(dit = {{54}, {112, 58}, {123, 102, 110}}; 

        n = Length[dit] + 1; 

        lowdiag = 
          Join[{Table[0, {n}]}, 
            Table[Join[dit[\([i]\)], Table[0, {n - Length[dit[\([i]\)]]}]], {
                i, n - 1}]]; 
d = lowdiag + Transpose[lowdiag]\)\(j = IdentityMatrix[n] - IdentityMatrix[n] + 1; 

      h = IdentityMatrix[n] - j/n; 
a = \(-d\)*d/2. ; 
b = h . a . h; 

      lambda = Eigenvalues[b]\  // \ Chop\)p = Count[lambda, x_ /; x > 0]
  charvect = Take[Eigenvectors[b], p];
points = Take[Sqrt[lambda], p]*charvect  //  Chop;
   Transpose[points] . points   pot = Transpose[points]Table[
            Apply[Plus, (pot[[r]] - pot[[s]])^2], {r, n}, {s, n}]
              //  Sqrt  //  Chop

ListPlot[pot, 
        Prolog -> PointSize[ .02]];

Timing[pot3 = makeconfiguration[dit3, 4]; ]dhat3 = newdistancediag[pot3]; stress[dit3, dhat3]

Markowitz portfolio

Optimal Portfolio

1,...,n     assets with random returns
R (nx1)  vector of expected returns
Σ (nxn)  covariance matrix of returns
X (nx1)  vector of weights of n assets  in portfolio
X.Σ.X    squared risk of the portfolio

(∑_(i = 1)^n) (Subscript[x, i]) = 1

I. Find X such that the expected return of the portfolio is prescribed, X.R=(Subscript[r, 0]), and the risk is minimal, short sales allowed, i.e., no restriction on signs of X.

II. The same objective, but short sales are not allowed, i.e., X must be  >= 0.

Eight major Czech companies (measured by market capitalization) with  the expected daily returns

R

ses = Inverse[
            Sqrt[DiagonalMatrix[
                Map[covariance[[#, #]]&, Range[Length[covariance]]]]]]; 
        
  
ses   .  covariance .  ses[[1]]

\(TraditionalForm
    \`\( (*\ \ Short\ sales\ allowed\ *) 
      \(portfolio[mean_, covariance_, δ_] := 
      Module[{x1, x2, inv, ones = Table[1, {Length[mean]}]}, 
        
inv = Inverse[covariance]; 
        
x1 = inv . mean/\((ones . inv . mean)\); 
        
x2 = inv . ones/\((ones . inv . ones)\); 
        
δ*x1 + \((1 - δ)\)*x2 // Simplify]
        
\)\)\)

optportfolio = portfolio[mean, covariance, δ]
  
(Subscript[r, 0]) = 0.00498274; 
optdelta = Solve[optportfolio .  R == 0.00498274, δ]
optportfolio . R // Simplify
s1 = Flatten[optportfolio /. optdelta]


   
          v1 = opt1 . Σ . opt1

?MultiplierMethod

nn = Length[mean]
xar = ToExpression[
          Map[StringJoin[{"x"}, #]&, Table[ToString[i], {i, nn}]]]
ones = Table[1, {nn}]; 
  
f = (xar . covariance . xar/2 // Simplify); 
  
g = {xar . ones - 1, xar . mean - 0.00498274}; 
  
h = {-xar} // Flatten
  
          
xar0 = Table[1. /nn, {nn}];

\(TraditionalForm
    \`\(sol = MultiplierMethod[f, g, h, xar, xar0]
\)\)

s2 = xar /. sol[[2]]
s2 . R

s1 . R

Options - Finance Package

Needs["Finance`Master`"];

?Finance`Options`*

This is a put option for an asset named XYZ that expires on December 20, 1992 with an exercise (or striking) price of 40.

xyzft = Option[{"Put", XYZ, 40, {12, 20, 1992}}]

Price[XYZ] ^= 40;
Volatility[XYZ] ^= 3/10;
(* Riskless interest rate *) r = 12/100;
(* The settlement date to be August 20, 1992 *)
settlement = {8, 20, 1992};
val1 = Value[xyzft, settlement, r]
N[val1]

?ImpliedVolatility

ImpliedVolatility[xyzft, settlement, r, val1] // 
      Chop

Plot[
        ImpliedVolatility[xyzft, settlement, r, v], {v, 1, 
          3}];