MVA II (HU Berlin, ERASMUS, July 3 - 12, 2008)

This Mathematica notebook will be displayed on the lecturer's web page.

http : // www.karlin.mff.cuni.cz/~hurt/

Mathematica Basics

References Wolfram, Hurt

Here you can look how to read Mathematica notebooks without Mathematica.

www.wolfram.com

Also a very fruitful source of information not  restricted to the Mathematica system only and providing a variety of knowledge from assorted fields:

http : // demonstrations.wolfram.com/

General notes:

1. After typing a command, start every calculation with  <Shift> + <Enter> or just <Enter> on the numeric part of the keyboard. <Enter> will just return to the next row in the input cell.

2. Built-in functions always begin with a capital. The system is case sensitive. Arguments of function are in [ ] brackets, algebraic priorities are emphasized by ( ), and { } are used for iterators ( {x,0,π} or {i,1,n}) or lists ( = vectors, matrices, trees, tables, &c.)

3. the cell (marked by the blue bracket on the right are basically input and output. Every input cell is introduced as In[34]:= while the corresponding output cell (after <Shift> + <Enter>) is Out[34]= .

4. If you do not need to see the output cell, it is a good idea to end the command with semicolon. The calculation is done but the output is not displayed. Appropriate for obvious results of the input or for large amounts of data in the output. 

5. Put the cursor on or in the command you'd like to explain and press F1.

4 pillars of Mathematica

Numerics

Simple calculator:

2+2

π is the symbol:

Sin[(π/2)]

Numerical value of π with approx. 50 digits:

N[π,50]

Graphics

Simple:

Plot[ Sin[(x^2)],{x,0,2π}]

More functions in one plot :

Plot[{(x^2) Sin[(x^2)],x Sin[2(x^2)], Sin[(x^2)]},{x,0,π},PlotRange->All]

Three-dimensional graphics:

Plot3D[Sin[(x^2)+y],{x,0,π},{y,0,π}]

ContourPlot[Cos[x]+Cos[y],{x,0,4Pi},{y,0,4Pi}]

DensityPlot[Sin[x]Sin[y],{x,-4,4},{y,-3,3},ColorFunction->"SunsetColors"]

Symbolics

e1=(((x+y))^11)

e2=Expand[e1]

Factor[e2]

The polynomial in the following equation may be factored, so explicit solution is achieved.

Solve[(x^5)-1==0,x]

Generally, by Solve you cannot obtain the solution if the polynomial is of order higher than 4.

Solve[(x^5)+(x^2)+1==0,x]

Nevertheless, you can obtain all the roots use the numeric equivalent to Solve:

sol5=NSolve[(x^5)+(x^2)+1==0,x]

For reasons that come over the scope of this lecture, we do not explain a perhaps rather strange form of the output. It is a general feature of Mathematica,  that the results are in the form of substitutions rules. Since Mathematica is in some way a rewriting language, the above form of the output has the advantage that it does not change the value of any variable involved. Nevertheless, there is a simple method to get the solution in the form of a list (roughly speaking, the symbol  "/." means substitute in accordance with the substitution rules.

roots5=x/.sol5

Further, we can extract the first, the second, and the last element of the above list by the following construction (we will put them into the list).

{First[roots5],roots5[[2]],Last[roots5]}

By Solve or NSolve we are able to solve the equations of the algebraic form. We will show it on the calculation of the internal rate of return (IRR). Suppose that c={(Subscript[c, 0]),(Subscript[c, 1]),...,(Subscript[c, n])} is a cash flow in equally spaced time intervals and i is the interest rate. Then the present value is

PV(c,i)=(∑_(j=0)^n)((Subscript[c, j])/(((1+i))^j))

and IRR is defined as a positive solution (with respect to i) to the equation

PV (c, i) = 0

First define the present value:

cf={-1050,100,100,100,100,1100};
presentvalue[i_]:=cf.(((1+i))^(-(Range[6]-1)))
presentvalue[i]

Next we solve the equation. We get also the complex roots which have no economic meaning.

irr=NSolve[presentvalue[i]==0,i]

The following construction returns only the positive roots

Cases[i/.irr,_Real]

Consider the following function and try to find its roots.

Plot[(E^x)-Cos[x],{x,-3π,π}]

The equation

eq1=(E^x)-Cos[x]==0

cannot be solved with Solve or NSolve. We get the warnings.

Solve[eq1]
NSolve[eq1]

In this case we can use another Mathematica built-in function, FindRoot. In this case we must give an initial value. We must be careful since for different initial values we obtain different solutions in this specific case.

FindRoot[(E^x)-Cos[x],{x,-2}]

FindRoot[(E^x)-Cos[x],{x,1}]

FindRoot[(E^x)-Cos[x],{x,-5}]

Examples of integration:

Integrate[Cos[x],x]

Integrate[(E^(-(x^2)/2)),{x,0,∞}]

If the indefinite integral may be expressed in terms of known special functions, Mathematica does it.

Integrate[Sin[(x^2)],x]

Integrate[Sin[(x^(3/2))],{x,0,π}]

Integrate[Sin[(x^3)],{x,0,π}]

This is the case where Mathematica does not know the answer and thus returns the input formula:

Integrate[Sin[Sin[(x^2)]],{x,0,π}]

Still we can get the numeric answer :

NIntegrate[Sin[Sin[(x^2)]],{x,0,π}]

Programming

A function definition that will take a single argument (the underscore _ stands for the formal argument in the definition):

hilb[n_]:=Table[(1/(i+j-1)),{j,1,n},{i,1,n}]

hilb[5]//MatrixForm

Multidimensional Scaling (MDS)

Remark

Since the built-in Mathematica functions and also the functions of AddOns packages begin always with capital letter, it is a recommended practise to define your own functions with the leading small letter. To distinguish between matrices and other quantities, I almost always use doubled symbol for a matrix, so that instead of  X, I write xx.

Demonstration that H.A.H  might be not p. s. d.

Recall that D is the distance matrix. It means that D is nonnegative symmetric with zeros on the diagonal. Matrix H is the centering matrix, see the very readable definition below. The reason for H.A.H might  not be positive semidivine is that A=-(1/2)D*D is not a matrix product but the Hadamard product (elementwise product of two matrices of the same dimensions)!!! In case of matrix product instead, matrix B must be positive semidefinite!

The following procedures are of general importance and should be carefully investigated. Function rd makes use of the idea that for arbitrary matrix X, the matrix X.Xᵀ is a positive semidefinite matrix. This is very useful piece of observation in situations if we need any positive semidefinite matrix for some testing purposes. It is very difficult to construct such a matrix by hand.

rd[n_]:=With[{xx=RandomReal[{0,1},{n,n}]},xx.Transpose[xx]-DiagonalMatrix[Diagonal[xx.Transpose[xx]]]]
centeringmatrix[n_]:=IdentityMatrix[n]-(1/n)Table[1,{n},{n}]
bb[n_]:=With[{dd=rd[n],hh=centeringmatrix[n]},-(1/2)hh.(dd*dd).hh]

Repeat the following command till you will see obviously negative eigenvalue.

bb[5]//Eigenvalues

Note that since the centering matrix H is idempotent of rank n-1, B can not have the rank greater than this number. So at least one of the eigenvalues is zero.

Implementation of the method

The first argument of the input is the distance matrix and the second (optional) p is the dimension of the configuration space. If p is greater than the number of positive eigenvalues of B then the dimension is set to the number of positive eigenvalues of B. If missing, then the dimension is set to the number of positive eigenvalues of matrix B. We make also use of "DistanceMatrix" function from the following standard Add-On package.

Needs["HierarchicalClustering`"];
Clear[msd];
msd[dd_(* the distance matrix *),p___(* possibly the desired dimension of the configuration space *)]:=
Module[{n,jj,aa,hh,bb,eigvec,eigval,eigsys,pdim,config,s1,ddhat},
n=Length[dd];
jj=Table[1,{n},{n}];
aa=-(1/2)dd*dd//N;
hh=IdentityMatrix[n]-(1/n)jj;
bb=hh.aa.hh;
eigsys=Reverse[Sort[Transpose[Eigensystem[bb]]]];
eigval=Map[First,eigsys];
Print["\nNumber of points = ",n,"\nEigenvalues:\n",eigval];
eigvec=Map[Last,eigsys];
If[p=!={},pdim=p];
pdim=Min[pdim,Count[eigval//Chop,_?Positive]];
eigval=Take[eigval,pdim];
eigvec=Take[eigvec,pdim];
Print["\n\nThe configuration in ",pdim,"-dimensional space is \n",(config=Transpose[Sqrt[eigval]*eigvec])//TableForm];
ddhat=DistanceMatrix[config]//Sqrt;
Print["Stress 1 (SS of the differences between original and estimated) = ",s1=Total[((Flatten[dd-ddhat])^2)]];
Print["Stress 2 (The square root of SS of the differences between original and estimated divided by the SS of original) = ",(Sqrt[(s1/(Total[((Flatten[dd])^2)]))])];
config
]

Sets of Data

Example 1

Distances of 4 Czech cities

(({{0, Praha-Podebrady, Praha-Hradec, Praha-Jihlava}, {\[Placeholder], 0, Podebrady-Hradec, Podebrady-Jihlava}, {\[Placeholder], \[Placeholder], 0, Hradec-Jihlava}, {\[Placeholder], \[Placeholder], \[Placeholder], \[Placeholder]}}))

distancematrix1=(({{0, 54, 112, 123}, {54, 0, 58, 102}, {112, 58, 0, 110}, {123, 102, 110, 0}}));

msd[distancematrix1]

Here is the plot of the resulting configuration:

ListPlot[{{-51.417108124407065`,-38.40597969279097`},{-24.01183114173669`,8.127769480186084`},{4.244989655753748`,58.78317163732917`},{71.18394961038999`,-28.504961424724318`}},PlotStyle->PointSize[0.02]]

Example 2

Given some points in (Subscript[E, 5]) but the calculated distance matrix is based on (Subscript[L, 1]) norm (Manhattan distance). The configuration is real but the input distance matrix is not Euclidean.

Needs["HierarchicalClustering`"]
(distancematrix2=DistanceMatrix[{{0,1,9,3,7},{1,0,8,2,6},{9,8,0,6,2},{3,2,6,0,4},{7,6,2,4,0},{1,2,3,4,5},{11,12,13,14,15}},DistanceFunction->ManhattanDistance]//N)//MatrixForm

Here we give the results for the desired dimension varying from 2 to 5. In Mathematica, it suffices to repeat the procedure with the proper number of dimensions in the form of the Table.

Table[msd[distancematrix2,j],{j,2,5}]

7 selected stocks (2 Frankfurt, 5 NYSE)

Reading and Exporting the Data to Stocks7.CSV

Real data from Mathematica Integrated Data Sources, till approx. July 9, 2008. If you are not connected to the Internet you can get the most actual data or play with the older data just changing "Jan. 1, 2008" to "Jan. 1, 2001" , say.   If you are not connected to the Internet just skip this Section and continue with the next Section Importing Data from the External File Stocks7.CSV

(* Do not calculate if you are not connected to the Internet!!! *)
priceIBM=FinancialData["IBM","Jan. 1, 2008"];
priceGE=FinancialData["GE","Jan. 1, 2008"];
priceBMW=FinancialData["F:BMW","Jan. 1, 2008"];
priceGM=FinancialData["GM","Jan. 1, 2008"];
priceJNJ=FinancialData["JNJ","Jan. 1, 2008"];
priceVWS=FinancialData["F:VWS","Jan. 1, 2008"];
priceABB=FinancialData["NYSE:ABB","Jan. 1, 2008"];

(* Do not calculate if not connected !!! *)prices={priceIBM,priceGE,priceBMW,priceGM,priceJNJ,priceVWS,priceABB};

Exporting prices to file "Stocks7.CSV"

(* Do not calculate if not connected !!! *)Export["C:\\Dokumenty\\VPFPM\\Berlin2008\\Stocks7.CSV",prices]

Importing Data from the External File Stocks7.CSV

im=Import["C:\\Dokumenty\\VPFPM\\Berlin2008\\Stocks7.CSV"];

prices={priceIBM,priceGE,priceBMW,priceGM,priceJNJ,priceVWS,priceABB}=Map[ToExpression,im,∞];

Showing how you can get parts of the lists from Mathematica.

prices//Last//First

priceABB//First

priceABB//First//First

Map[Last,prices]

Plotting the Prices

Map[DateListPlot[#,Joined->True]&,prices]

Find a list of NYSE stocks whose ticker symbols begin with IB :

FinancialData["NYSE:AB*","Lookup"]

Find a list of Frankfurt stocks whose ticker symbols begin with BM :

FinancialData["F:BM*","Lookup"]

Consolidation of the Data

We use data from two Stock Exchanges and the trading dates do not coincide. To compare the prices and/or returns we must first select the common trading dates and then extract the prices correspondingly. Check separately the parts of the following rather complex commands. You must also look how pure functions (# &) work, see Help.

commondates=Apply[Intersection,Map[First[Transpose[#]]&,prices]];
commondates//Last
commondates//Length

Here are two constructions how to extract the date and price of the IBM stock of April 18, 2008.

Select[priceIBM,First[#]=={2008,4,18}&]
Select[priceIBM,Function[x,First[x]=={2008,4,18}]]

By this database construction we obtain the database containing the dates and prices that are quoted on both Exchanges on the same dates.

consolidateddata=Map[Function[y,Map[Select[y,Function[x,First[x]==#]]&,commondates]],prices];

Number of companies:

consolidateddata//Length

A shortened output to see the structure of the data.

Short[consolidateddata,10]

For a further need we only need prices. Here is the method how to get the prices of the first stock, i. e. IBM:

Map[Last,Flatten[consolidateddata//First,1]]

And similarly, by a little bit more complicated construction, we get the prices of all 7 companies. The command "Dimensions" shows that from now we deal with the matrix with 7 rows representing stocks and 129 columns representing trading days.

consolidatedprices=Map[Map[Last,Flatten[consolidateddata[[#]],1]]&,Range[Length[consolidateddata]]];
Dimensions[consolidatedprices]

The Means, Covariances, and Correlations (prices)

Mean[consolidatedprices//Transpose]

Covariance[consolidatedprices//Transpose]//MatrixForm

Correlation[consolidatedprices//Transpose]//MatrixForm

Testing Random Walk Hypothesis

Next time.

Returns and their Characteristics

(* Returns given in per cent *)returns=100*Map[(Rest[consolidatedprices[[#]]]-Most[consolidatedprices[[#]]])/Most[consolidatedprices[[#]]]&,Range[Length[consolidateddata]]];

Map[ListPlot[#,Joined->True]&,returns]

Expected returns estimated from the historical time series. The maximum expected returns is for IBM, position 1.

(expectedreturns=Mean[returns//Transpose])//TableForm
{Max[expectedreturns],Position[expectedreturns,Max[expectedreturns]]}

(covariancematrixofreturns=Covariance[returns//Transpose])//MatrixForm

Now look on the expected return and risk (= standard deviation), all in per cent. We see that the displayed results are in the contradiction with the assumptions of the efficient market. The most risky is GM even with the highest negative expected return, but also other some of the assets with negative returns are more risky than those with the positive return.

Expected returns and standard deviations

{expectedreturns,covariancematrixofreturns//Diagonal//Sqrt}//Transpose//TableForm

(correlationmatrixofreturns=Correlation[returns//Transpose])//MatrixForm

Maximum correlation is between IBM and GE returns:

With[{c=correlationmatrixofreturns},{Max[c-IdentityMatrix[Length[c]]],Position[c,Max[c-IdentityMatrix[Length[c]]]]}]

Eigenvalues of the covariance matrix

covariancematrixofreturns//Eigenvalues

Eigenvalues of the correlation matrix

correlationmatrixofreturns//Eigenvalues

Markowitz portfolio, short sales allowed

Lagrange function generally, μ prescribed expected return on portfolio

Clear[lagrange];
lagrange[xx_,rr_,vv_,λ1_,λ2_,μ_]:=
Module[{boldone},
boldone=Table[1,{Length[xx]}];
(1/2)xx.vv.xx+λ1(1-boldone.xx)+λ2(μ-rr.xx)]

Lagrange function for the problem of 7 stocks:

Clear[x,xx,λ1,λ2];
xx=Array[x,{Length[expectedreturns]}];
ones=Table[1,{Length[expectedreturns]}];
lagrange7=lagrange[xx,expectedreturns,covariancematrixofreturns,λ1,λ2,μ]//Simplify

The gradient and the derivatives with respect to λ1 and λ2

(dxx=D[lagrange7,{xx}])//TableForm
dλ1=D[lagrange7,λ1]
dλ2=D[lagrange7,λ2]

Construction of the system of N+2 equations (N being the number of assets in question), here 7 + 2 =9. Still μ is not specified.

equations={Thread[dxx==0],dλ1==0,dλ2==0}//Flatten

solution=Solve[equations,{xx,λ1,λ2}//Flatten]

Optimal portfolio for an arbitrary μ:

optimalportfolio=(xx/.solution)//First

Check that the total investment is actually equal 1 (i. e. (∑_(n=1)^N)(Subscript[x, n])=1):

Total[optimalportfolio]

Now some concrete examples. The following table shows (μ,(Subscript[r, p]),(Sqrt[var (Subscript[ρ, p])]),x):

Table[{"μ:",m," Achieved",expectedreturns.optimalportfolio/.μ->m,"\n","Risk: ",√(optimalportfolio.covariancematrixofreturns.optimalportfolio)/.μ->m,"Portfolio: \n\n\n",optimalportfolio/.μ->m,"\n\n"},{m,0.04,0.14,0.02}]//TableForm

Remind the expected returns:

expectedreturns//TableForm

We see that short sales appear in all minimum-variance portfolios. With short sales allowed, the optimal portfolio may have the expected return higher than the maximum expected return of any individual asset in question.

Markowitz portfolio, short sales not allowed

This is the case of quadratic programming. We make use of the built - in function NMinimize. Notation as in the previous subsubsection. The error message (in red) relates to the last case where the prescribed return is 0.14 and it exceeds the maximum of available returns. We also see that with the inequality in the wealth condition, the total health is not completely spent.

Table[{"Prescribed return: ",μ," Portfolio: ",xx/.(NMinimize[{xx.covariancematrixofreturns.xx,xx.ones≤1,xx.expectedreturns≥μ,Thread[xx≥0]}//Flatten,xx]//Last),"\n\n"},{μ,0.04,0.14,0.02}]//TableForm

So if you are not too ambitious, you even do not spend your whole wealth. Minimizing risk is the prevailing criterion. If you really want to spend all of your money and you are also willing to avoid the problems with the equality constraints, you may allow your investment to be in a very narrow interval, [0.9999999, 1], say:

Print["Portfolio: ", x689=xx/.(NMinimize[{xx.covariancematrixofreturns.xx,0.9999999≤ xx.ones≤ 1,xx.expectedreturns≥0.12,Thread[xx≥0]}//Flatten,xx]//Last),"with the sum of weights: ",
Total[x689]]

Compare with the above inequality constraint: the total wealth of 1 is not spent.

Print["Portfolio: ", x689=xx/.(NMinimize[{xx.covariancematrixofreturns.xx, xx.ones≤ 1,xx.expectedreturns≥0.12,Thread[xx≥0]}//Flatten,xx]//Last),"with the sum of weights: ",
Total[x689]]

There might be problems with the equality constraint. In some small examples we likely do not meet them.

Print["Portfolio with the prescribed return 0.1 is ",xx/.(NMinimize[{xx.covariancematrixofreturns.xx,xx.ones==1,xx.expectedreturns≥0.1,Thread[xx≥0]}//Flatten,xx]//Last)]

Print["Portfolio with the prescribed return 0.13 is ",xx/.(NMinimize[{xx.covariancematrixofreturns.xx,xx.ones==1,xx.expectedreturns≥0.13,Thread[xx≥0]}//Flatten,xx]//Last)]

Sharpe ratio, short sales  allowed

We use function NMaximize without any constraints, and then just normalize the solution:

Module[{solutionSharpeShort},solutionSharpeShort=xx/.(NMaximize[((expectedreturns.xx)/(Sqrt[xx.covariancematrixofreturns.xx])),xx]//Last);
Print["The optimal portfolio is\n",(solutionSharpeShort/(Total[solutionSharpeShort])),"\n and it really sums to one: ",(solutionSharpeShort/(Total[solutionSharpeShort]))//Total]]

Sharpe ratio, short sales not allowed

To exclude very small (from numerical point of view zeros) we use the built-in function Chop.

Module[{solutionSharpe},solutionSharpe=solutionSharpe=xx/.(NMaximize[{((expectedreturns.xx)/(Sqrt[xx.covariancematrixofreturns.xx])),Thread[xx≥0]}//Flatten,xx]//Last);
Print["The optimal portfolio is\n",(solutionSharpe/(Total[solutionSharpe]))//Chop,"\n and it really sums to one: ",(solutionSharpe/(Total[solutionSharpe]))//Total]]