Foundations of Computational Finance

Ramesh Adhikari

The Wolfram Language has numerous knowledge-based built-in functions to support financial computations. This article introduces many built-in and other financial functions that are based on concepts and models covered in undergraduate-level finance courses. Examples are taken from a wide range of finance areas. They emphasize importing and visualization of data from many sources, valuation, capital budgeting, analysis of stock returns, portfolio optimization and analysis of bonds and stock options. We hope that all the functions selected in this article are very useful for analyzing real-world financial data. All examples provide a unique set of tools for users to engage with real-world financial data and solve practical problems. The feature of automatic data retrieval from online sources and its analysis makes all results reproducible without any modifications in the code. We hope this feature will attract new users from the finance community.

1. Introduction

Finance is computational in nature and often involves the analysis and visualization of complex data, optimization, simulation and use of data for risk management. Without the proper use of technology, it is almost impossible to analyze these functions of modern finance. Moreover, the field of finance has become far more driven by data and technology since the 1990s, which has made large-scale data analysis the norm. Data-driven decision-making and predictive modeling are now the heart of every strategic financial decision. Since the publication of Varian [, ], Shaw [] and Stojanovic [], there have been many updates and new functions, but no new articles or books have been written to cover wide areas of computational finance. This article provides a comprehensive overview of functions related to finance and introduces many functions that are useful for real-world financial data analysis using Mathematica 12.

We have provided all the custom functions in the text so that users can make changes as they learn how to program in the Wolfram Language. Furthermore, we minimize the explanation of any financial concepts in this article as our focus is on introducing financial application of the Wolfram Language.

We begin by defining some symbols that are frequently used as input arguments of custom functions in the article. The article uses symbol or symbols, startDate, endDate and period as arguments in many functions defined in this article. Most of these symbols are used as input in the built-in function FinancialData. All these arguments must be specified in the format acceptable in the FinancialData function. We use symbol or symbols to represent a company’s or companies’ stock ticker symbol or symbols. It could be a string or a list of strings. The format startDate represents the start date of the sample period specified and endDate represents the last date of the analysis period. Both must be specified as date objects in any date format supported by DateObject. Similarly, period represents data frequency. It may include “Day”, “Week”, “Month” or “Year”. In the subsequent functions defined in this article, we will not describe them when they are used as arguments.

The article is organized into 13 sections:

1. this Introduction

. importing and visualizing data from different sources

. capital budgeting and business valuation

. functions for the analysis of security returns

. rolling-window performance analysis

. financial application of optimization

. decomposing the risk of a portfolio into its components

. importing factor data and running factor models

. computing different types of portfolio performance measures

. technical analysis of stock prices

. bond analysis

. analyzing derivative products

. concluding remarks

2. Accessing Financial Data and Basic Data Visualization

The most commonly used built-in functions for retrieving company-specific financial data are Company, CompanyData and Financial. For example,

CompanyData[EntityValue[Interpreter["Company"]["FB"],"Company"],"Dataset"]

imports Facebook's financial statement data and

Interpreter["Financial"]["FB"]["Dataset"]

imports Facebook's price-related data.

Similarly, the function FinancialData can be used to get data about stocks and other financial instruments. DateListPlot or CandlestickChart can be used to chart prices against time. TradingChart or InteractiveTradingChart can be used to make interactive plots with additional features of adding different technical indicators. Other functions such as KagiChart, RenkoChart, PointFigureChart or LineBreakChart can also be used to visualize financial data. In the remaining part of this section, we are going to show you how to import data from different sources and visualize it.

plotCompanyData[symbol_String,startDate_,endDate_,properties:{__String}]:=
DateListPlot[
EntityValue[
EntityValue[Interpreter["Financial"][symbol],"Company"],Dated[properties,Interval[{DateObject[startDate],DateObject[endDate]}]]
],
PlotLabel->symbol,PlotLegends->Placed[properties,Bottom],PlotRange->All
]

We download Apple’s return on assets (ROA), return on equity (ROE) and revenue growth over the period January 1, 2001, to January 1, 2019, and plot them.

plotCompanyData["AAPL",{2000,1,1},{2019,1,1},{"ROA","ROE","RevenueGrowth"}]

Similarly, we define the plotCompaniesData function to compare any specified property of different companies. The function takes a list of stock symbols, beginning period, end period and a property to consider as its arguments.

plotCompaniesData[symbols:{__String},startDate_,endDate_,property_String]:=
DateListPlot[
EntityValue[
EntityValue[Interpreter["Financial"][#],"Company"]&/@symbols,
Dated[property,Interval[{DateObject[startDate],DateObject[endDate]}]]
],
PlotLabel->property,PlotLegends->Placed[symbols,Bottom],PlotRange->Full,PlotRange->All
]

We plot the revenue growth of Apple, Facebook, Walmart and Bank of America over the period January 1, 2000, to January 1, 2019.

plotCompaniesData[{"AAPL","FB","WMT","BAC"},{2000,1,1},{2019,1,1},"RevenueGrowth"]

Second, we import and visualize data from the Federal Reserve Bank of St. Louis, as it is one of the most important data sources when it comes to economic data. The built-in function ServiceExecute can be used to request the data from the Federal Reserve Economic Data API. Its argument structure is:

ServiceExecute["FederalReserveEconomicData","SeriesData","ID"->{ids}]

where ids is a series ID or a list of IDs. It returns a time series containing data for the specified series. It is often of interest to plot the economic time series with the recession dates.

The plotFedData function downloads and plots the selected series along with the shaded recession period. The function takes series ID, start date, end date and title as inputs and returns a graph. It uses recession indicators based on USREC (US recession) data from the National Bureau of Economic Research (NBER) for the United States from the period following the peak through the trough to indicate the recession period. The Federal Reserve Bank of St. Louis may require the API key to download its data. The API key can be obtained freely by creating a user account at https://fred.stlouisfed.org (click "my account" and follow the instructions).

plotFedData[property_String,startDate_,endDate_,title_String]:=Module[
{conn,recessionIndicators,recessionDates,data},
conn=ServiceConnect["FederalReserveEconomicData"];
recessionIndicators=TimeSeriesWindow[conn["SeriesData","ID"->{"USREC"}],{startDate,endDate}];
recessionDates=Join[recessionIndicators["Dates"][[#]]&/@Position[Prepend[Differences[recessionIndicators["Values"]],0],1.],recessionIndicators["Dates"][[#]]&/@Position[Append[Differences[recessionIndicators["Values"]],0],-1.],2];
data=TimeSeriesWindow[conn["SeriesData","ID"->{property}],{startDate,endDate}];
Show[
DateListPlot[
Thread[{#,Min[data["Values"]]}]&/@recessionDates,
Joined->True,Filling->Top,FillingStyle->Opacity[0.3,Gray],PlotStyle->None,PlotRange->{{First[data["Dates"]],Last[data["Dates"]]},{Min[data["Values"]],Max[data["Values"]]}},PlotTheme->"Detailed",PlotLabel->title],
DateListPlot[data,Joined->True,PlotRange->All]
]
]

Now we can download any series and plot it. For example, we download and plot the Leading Index for the United States (USSLIND) over the period January 30, 1990, to January 30, 2019. Please use the API key 207071a5f2e90e7816259d3c32c1ab81 if needed. The shaded regions indicate recession periods.

plotFedData["USSLIND",{1990,01,30},{2019,01,30},"Leading Index for United States"]

getQuandlData[dataID_String,startDate_,endDate_,APIKey_String]:=Module[
{data},
data=Import[
StringJoin[
"https://www.quandl.com/api/v3/datasets/",
dataID,
"/data.csv?start_date=",
startDate,
"&end_date=",
endDate,
"&order=asc&api_key=",
APIKey]
];
TimeSeries[data[[2;;,2]],{ToExpression[StringSplit[data[[2;;,1]],"-"]]}]
]

We download and plot the historical real S&P 500 prices by month (MULTPL/SP500_REAL_PRICE _MONTH) over the period March 31, 1975, to May 30, 2019.

DateListPlot[
getQuandlData["MULTPL/SP500_REAL_PRICE_MONTH","1975-03-31","2019-05-30","G_ydzsMB4zEWcWz3xNs4"],
GridLines->Automatic,
PlotLabel->"Monthly S&P 500 Real Price"
]

Finally, we show how to create a dataset. The built-in function Dataset is very useful for organizing large or small sets of data.

The function getCompanyData can be used to get a company's fundamental data. After the data is stored, we can organize and analyze the data.

getCompanyData[symbol_String,startDate_,endDate_,properties:{__String}]:=Module[
{data,mergedData,dates},
data=EntityValue[
EntityValue[Interpreter["Financial"][symbol],"Company"],Dated[properties,Interval[{DateObject[startDate],DateObject[endDate]}]]
];
mergedData=TimeSeriesThread[#&,data,ResamplingMethod->"NaN"];
dates=DateObject[#]&/@mergedData["Dates"][[All,1,{1,2,3}]];
Dataset@Map[
AssociationThread[Flatten[{"Date",properties}]->#]&,Transpose@Insert[Transpose[QuantityMagnitude[mergedData]["Values"]],dates,1]
]
]

For example, we download return on assets ("ROA"), return on equity ("ROE") and revenue growth ("RevenueGrowth") for Apple Inc. (AAPL) over the period January 1, 2000, to January 1, 2019, and make a dataset. After the dataset is constructed, we can pull data and do further analysis using a rich set of built-in Wolfram knowledge.

getCompanyData["AAPL",{2000,1,1},{2019,1,1},{"ROA","ROE","RevenueGrowth"}]

3. Common Financial Decision Making, Capital Budgeting and Business Valuation

Basic concepts used in common financial decision making are important for learning and understanding the finance discipline. Many functions such as EffectiveInterest, CashFlow, Annuity, AnnuityDue and TimeValue are directly related to finance. Other functions such as Solve and FindRoot can be used to find one of the unknowns when relevant information is given. All these functions are useful for solving time value of money, capital begetting and business valuation problems. The Mathematica documentation provides numerous examples of how to use these functions. In this section, we are going to focus on a few examples concerning loan amortization, capital budgeting and business valuation.

A loan amortization table is often used to visualize periodic payments of the loan, loan balance and payment breakdown into principle payment and interest payment. The function amortizationTable returns an amortization table given its input arguments. It takes four arguments:

pv: current value of loan amount

years: loan term in years

APR: annual percentage interest rate

frequency: frequency of loan payment per year: 12 for monthly payment, 1 for annual payment, and so on

fv: (an optional argument) future value of the loan amount; if no value for fv is provided, the future value of the loan is assumed to be zero

amortizationTable[pv_,years_,APR_,frequency_,fv_:0]:=Module[
{periods,rate,payment,balance,interest,principal},
periods=years frequency;
rate=APR/frequency;
payment=(pv+((pv+fv)/((((1+rate))^periods)-1)))rate;
balance=TimeValue[Annuity[payment,periods],rate,0];
Text@Grid[
Prepend[
Round@Table[
interest=balance rate;
principal=payment-balance rate;
balance=balance-principal;
{month,balance+principal,payment,interest,principal,balance},
{month,periods}
],
{"Period","Beginning Balance","Payment","Interest Payment","Principal Payment","Ending Balance"}
],Dividers->All]]

Using this function, we compute an amortization table for a loan of $40,000 with 1-year loan term, 5% APR paid monthly.

amortizationTable[40000,1,0.05,12]

The most commonly used decision tools in capital budgeting are net present value (NPV), internal rate of return (IRR), modified internal rate of return (MIRR) and profitability index (PI). These are defined in terms of the cash flows (Subscript[CF, i]), i=0,1,2,…,n, the discount rate r and the reinvestment rate rr by:

NPV=((Subscript[CF, 1])/((1+r)))+((Subscript[CF, 2])/(((1+r))^2))+…+((Subscript[CF, n])/(((1+r))^n))+(Subscript[CF, 0]),

((Subscript[CF, 1])/((1+IRR)))+((Subscript[CF, 2])/(((1+IRR))^2))+…+((Subscript[CF, n])/(((1+IRR))^n))+(Subscript[CF, 0])=0,

((((Subscript[CF, 1])(1+rr))^(n-1))+(((Subscript[CF, 2])(1+rr))^(n-2))+…+(Subscript[CF, n-1])(1+rr)+(Subscript[CF, n]))/(((1+MIRR))^n)=(Subscript[CF, 0]),

PI=((((Subscript[CF, 1])/((1+r)))+((Subscript[CF, 2])/(((1+r))^2))+…+((Subscript[CF, n])/(((1+r))^n)))/(Subscript[CF, 0])).

We use the built-in Mathematica functions TimeValue, FindRoot and Cashflow in the capitalBudgetingTools function to compute these measures. It takes cash flows CF (a list), discount rate r and reinvestment rate rr as its arguments.

capitalBudgetingTools[CF_List,r_,rr_]:=
{
"NPV"->TimeValue[Cashflow[CF],r],
First@FindRoot[TimeValue[Cashflow[CF],IRR]==0,{IRR,0.005}],
First@FindRoot[TimeValue[TimeValue[Cashflow[Flatten[{0,CF[[2;;]]}]],rr,Length[Rest[CF]]],MIRR,-Length[Rest[CF]]]==-First[CF],{MIRR,0.005}],
"PI"->TimeValue[Cashflow[Flatten[{0,CF[[2;;]]}]],r]/Abs[First[CF]]
}

We illustrate the use of the capitalBudgetingTools function with an example. Say a project requires a $50,000 initial investment and is expected to produce, after tax, a cash flow of $15,000, $8,000, $10,000, $12,000, $14,000 and $16,000 over the next six years. The discount rate is 10% and the reinvestment rate is 11%. We compute the project’s NPV, IRR, MIRR and PI.

capitalBudgetingTools[{-50000,15000,8000,10000,12000,14000,16000},0.1,0.11]

One of the most widely used business valuation models is the discounted cash flow model, in which the value of any asset is obtained by discounting the expected cash flows on that asset at a rate that reflects its riskiness. In its most general form, the value of a company is the present value of the expected free cash flows the company can generate in perpetuity. Because we cannot estimate free cash flows (FCF) in perpetuity, we generally allow for a period where FCF can grow at extraordinary rates, but we allow for closure in the model by assuming that the growth rate will decline to a stable rate that can be sustained forever at some point in the future. If we assume that the discount rate is the weighted average cost of capital (WACC), FCF grows at the rate of g per year and that the last year’s free cash flow is (Subscript[FCF, 0])>0, then the value of the firm can be defined as

(Subscript[Value, 0])=(∑_(t=1)^∞)(((Subscript[FCF, 0]) (1+g) )/(((1+WACC))^t)).

If we assume that FCF grows at the rate (Subscript[g, 1]) for the next (Subscript[n, 1]) years and at the rate (Subscript[g, 2]) thereafter, then the value of the firm can be written as

(Subscript[Value, 0])=(∑_(t=1)^(Subscript[n, 1]))(((Subscript[FCF, 0]) (((1+(Subscript[g, 1])) )^t))/(((1+WACC))^t))+(1/(((1+WACC))^(Subscript[n, 1])))×(((Subscript[FCF, (Subscript[n, 1])]) (1+(Subscript[g, 2])) )/((WACC-(Subscript[g, 2])))).

We implement formula () with twoStageDCFModel that takes five arguments:

1. FCF, last year's free cash flows

2. g1, the annual growth rate of free cash flows in the first growth period

3. n1, the number of years in the first growth period

4. g2, the stable growth rate

5. wacc, the weighted average cost of capital

twoStageDCFModel[FCF_,g1_,n1_,g2_,wacc_]:=(∑_(t=1)^n1)(((FCF×(((1+g1))^t))/(((1+wacc))^t)))+((FCF×(((1+g1))^n1)×(1+g2))/((((1+wacc))^n1)×(wacc-g2)))

Suppose that a company has $100 as its past year’s CFC, its FCF is expected to grow at 5% for the next five years and 2.5% thereafter and that its WACC is 9.5%. We compute its value.

twoStageDCFModel[100,0.05,5,0.025,0.095]

The most practical approach is to assume that a company goes through three phases: growth, transition and maturity. First, a company’s growth increases. In the transition phase, the growth rate decreases. In the mature phase, a company grows at the same rate as that of the overall economy.

Assume that the FCF is positive and grew at the rate (Subscript[g, 0]) last year. Assume further that it grows at the higher rate (Subscript[g, 1]) each year for the next (Subscript[n, 1]) years and that after (Subscript[n, 1]) years, it declines at the rate (Subscript[g, 2]) each year for the next (Subscript[n, 2]) years. Finally, assume the company grows at the stable positive rate (Subscript[g, 3]) per year after (Subscript[n, 1])+(Subscript[n, 2]) years and that the cost of capital is represented by WACC. Then the value of the firm can be written as

(Subscript[Value, 0])=(∑_(t=1)^(Subscript[n, 1]))(((Subscript[FCF, t-1]) (((1+(Subscript[g, 0])+t × (Subscript[g, 1])) )^t))/(((1+WACC))^t))+(∑_(t=(Subscript[n, 1])+1)^((Subscript[n, 1])+(Subscript[n, 2])))(((Subscript[FCF, t-1]) (((1+(Subscript[g, 0])+(Subscript[n, 1])×(Subscript[g, 1])+(t-(Subscript[n, 1])) × (Subscript[g, 2])) )^((t-(Subscript[n, 1])))))/(((1+WACC))^t))+(1/(((1+WACC))^((Subscript[n, 1])+(Subscript[n, 2]))))×(((Subscript[FCF, (Subscript[n, 1])+(Subscript[n, 2])]) (1+(Subscript[g, 3])) )/((WACC-(Subscript[g, 3]))))

We implement formula () with threeStageDCFModel that takes eight arguments:

1. FCF, last year's free cash flows

2. g0, last year's FCF growth rate

3. g1, the incremental growth rate in the high growth period

4. n1, the number of years in the high growth period

5. g2, the declining growth rate in the translational growth period

6. n2, the number of years in the translational growth period

7. g3, the stable growth rate in the maturity growth period

8. wacc, the weighted average cost of capital

threeStageDCFModel[FCF_,g0_,g1_,n1_,g2_,n2_,g3_,wacc_]:=N@Module[
{growthrates,cashflows},
growthrates=1+Flatten[{g0+Range[n1] g1,Last[g0+Range[n1] g1]+Range[n2] g2}];
cashflows=FCF FoldList[Times,growthrates];
TimeValue[Cashflow[Flatten[{0,cashflows}]],wacc,0]+Last[cashflows] (1+g3)/(wacc-g3)/(1+wacc)^(n1+n2)
]

To apply the function, consider a company in the early stage of its life cycle, assuming that the company experienced 10 percent growth in the past year. The company is expected to grow by 8% more each year for the next 7 years and its growth will start to decline by 5% each year after the seventh year for 5 years. After 12 years, the company is expected to grow at the same rate as that of the overall economy, which is 2.5% per year. Suppose the past year’s CFC was $100 million and the weighted average cost of capital is 9.5%. We compute the value of the company.

threeStageDCFModel[100,0.1,0.08,7,-0.05,7,0.025,0.095]

When using threeStageDCFModel, we assume that the growth rate in the stable phase is never negative. Therefore, do not assume the declining growth rate to be too high. Otherwise, the growth rate of FCF in the maturity phase may be negative.

4. Analysis of Stock Returns

There are various methods in analyzing historical stock prices and returns data. We can use different kinds of charts and graphs as well as descriptive statistics. Similarly, it is also important to understand whether the stock returns distribution is normal. In the next two subsections, we explain some of the most common charts and descriptive measures.

4.1 Individual Stock Return Analysis

Commonly used charts for historical performance analysis are time series plot of prices, normalized prices (historical prices divided by the price at the beginning), continuous draw-downs (cumulative continuous returns) and cumulative returns. The function performanceOfOneSymbol takes four arguments as defined in Section  and returns four different plots: historical prices, the normalized price, continuous draw-downs and cumulative returns.

performanceOfOneSymbol[symbol_String,startDate_,endDate_,period_]:=Module[
{data,prices,dates,returns,timeseries,drawdowns,indices},data=FinancialData[symbol,"Close",{startDate,endDate,period}];
prices=QuantityMagnitude[data["Values"]];
dates=data["Dates"][[All,1,{1,2,3}]];
returns=100 Differences[prices]/Most[prices];
timeseries=TimeSeries[returns,{Rest[dates]}];
drawdowns=TimeSeries[Flatten[Accumulate/@SplitBy[returns,Sign]]/._?Positive->0,{dates[[2;;]]}];
indices=TimeSeries[prices/Table[First[prices],Length[prices]],{dates}];
{DateListPlot[data,PlotLabel->"Prices over Time",PlotTheme->"Detailed"],DateListPlot[indices,PlotLabel->"Normalized Price over Time",PlotTheme->"Detailed"],DateListPlot[Accumulate[timeseries],PlotLabel->"Cumulative Returns over Time",PlotTheme->"Detailed"],DateListPlot[drawdowns,PlotRange->All,Filling->Axis,PlotTheme->"Detailed",PlotLabel->"Continuous Draw-down"]}]

We can apply performanceOfOneSymbol to any symbol and period. For example, we plot the historical stock prices and returns of Walmart Inc. (WMT) over the period October 10, 2018, to June 7, 2019.

performanceOfOneSymbol["WMT",{2018,10,10},{2019,6,7},"Day"]

A histogram and an empirical plot of kernel density estimates are often used to describe the general shape of the data. The densityPlot function takes four arguments as defined in Section  and returns density and histogram plots of returns.

densityPlot[symbol_String,startDate_,endDate_,period_]:=Module[
{prices,returns,x},
prices=FinancialData[symbol,"Close",{startDate,endDate,period}];
returns=Differences[Log[QuantityMagnitude[prices]["Values"]]];
Show[
Histogram[returns,Automatic,"PDF",PlotLabel->"Daily Returns",PlotTheme->"Detailed"],
Plot[{
PDF[SmoothKernelDistribution[returns],x],PDF[NormalDistribution[Mean[returns],StandardDeviation[returns]],x]},{x,Min[returns],Max[returns]},PlotRange->All,PlotLegends->Placed[{"Empirical Density","Normal Density"},Below],PlotStyle->{Dashed,Thick}]
]
]

Using the densityPlot function, we download the daily closing price of the S&P 500 index over the period October 1, 2000, to October 1, 2019, and plot the histogram, the empirical density function and the density function of a normal distribution with the same mean and variance.

densityPlot["^SPX",{2000,10,1},{2019,10,1},"Day"]

Many built-in functions can be used to compute different descriptive statistics. These descriptive statistics describe properties of distributions, such as location, dispersion and shape. The most common measures are computed by performanceStatsOneSymbol:

• holding period return

• average return

• geometric mean return

• cumulative returns

• standard deviation

• minimum return

• maximum return

• skewness

• kurtosis

• historical value at risk

• historical conditional value at risk

performanceStatsOneSymbol[symbol_String,startDate_,endDate_,period_]:=Module[
{data,prices,dates,returns,var,cvar},
data=FinancialData[symbol,"Close",{startDate,endDate,period}];
prices=QuantityMagnitude[data["Values"]];
dates=data["Dates"][[All,1,{1,2,3}]];
returns=Differences[prices]/Most[prices];
var=Sort[returns][[Floor[Length[returns] 0.05]]];
cvar=Mean[Select[returns,#≤var&]];
Apply[#1->#2&,
Transpose[{
{
"Symbol",
"Start Date",
"End Date",
"Holding Period Return",
"Arithmetic Mean",
"Geometric Mean",
"Cumulative Return",
"Stdandard Deviation",
"Minimum Return",
"Maximum Return",
"Skewness",
"Kurtosis",
"Value at Risk (VaR)",
"Conditional VaR (CVaR)"
},
{
symbol,
FromDigits[First[dates][[1;;2]]],
FromDigits[Last[dates][[1;;2]]],
100 (Last[prices]/First[prices]-1),
100 Mean[returns],
100 (GeometricMean[1+returns]-1),
100 Total[returns],
100 StandardDeviation[returns],
100 Min[returns],
100 Max[returns],
Skewness[returns],
Kurtosis[returns],
100 var,
100 cvar
}
}],{1}]
]

For example, we compute the descriptive statistics for Walmart Inc. (WMT) using monthly returns over the period October 10, 2010, to June 6, 2019.

performanceStatsOneSymbol["WMT",{2010,10,10},{2019,6,6},"Month"]//Column//Text

It is also informative to examine the historical performance of an individual stock. In such an analysis, we calculate monthly statistics using daily returns and report them on a monthly basis. We define the monthlyOneStockStatistics function to download historical stock prices, compute desired statistics and return a dataset. The function takes four arguments: stock ticker symbol, start date, end date and a statistical function such as Mean, StandardDeviation or Total, and returns a dataset. For the function to work, it requires more than two years of data.

monthlyOneStockStatitics[symbol_String,startData_,endDate_,statFunction_]:=Module[
{data,prices,returns,dates,years,months,timeseries,results,firstRow,annualReturns,midRows,lastRow,heads},
data=FinancialData[symbol,{startData,endDate}];
prices=QuantityMagnitude[data["Values"]];
returns=100 Differences[prices]/Most[prices];
dates=data["Dates"][[All,1,{1,2,3}]];
years=DeleteDuplicates[DateString[#,"Year"]&/@dates[[2;;]]];
months=dates[[2;;]][[All,2]];
timeseries=Transpose@{DateString[#,"Year"]<>DateString[#,"Month"]&/@dates[[2;;]],returns};
results=statFunction/@GatherBy[timeseries,First];
firstRow=Append[
Join[
Table["NA",First[months]-1],
results[[1;;13-First[months],2]]
],
statFunction[results[[1;;13-First[months],2]]]
];
annualReturns=statFunction/@Partition[results[[(13-First[months]+1);;-(Last[months]+1),2]],12];
midRows=Transpose@Append[
Transpose[Partition[results[[(13-First[months]+1);;-(Last[months]+1),2]],12]],
annualReturns
];
lastRow=Append[
Join[
results[[-Last[months];;,2]],
Table["NA",12-Last[months]]
],
statFunction[results[[-Last[months];;,2]]]
];
heads={"Year","January","February","March","April","May","June","July","August","September","Octber","November","December","Annual"};
Dataset[
AssociationThread[heads->#]&/@
Transpose[Insert[Transpose@Append[Prepend[midRows,firstRow],lastRow],years,1]]
]
]

For example, we compute the monthly cumulative returns for Walmart Inc. (WMT) using daily returns over the period October 1, 2010, to June 30, 2019.

monthlystatistics=monthlyOneStockStatitics["WMT",{2010,10,1},{2019,06,30},Total];

Once we compute the statistics, we can take specific columns by specifying their names.

monthlystatistics[All,{"Year","January","February","March","Annual"}]

We are often interested in knowing whether returns data follows a normal distribution because understanding whether stock returns are normal or not is very important in investment management. One way to check whether returns are normally distributed or not is to compare the empirical quantiles of the data with normal distribution. The QuantilePlot function can be used to produce quantile-quantile plots. Many other built-in functions can help to assess whether returns are normally distributed.

The DistributionFitTest function can be used to test whether data is normally distributed and can also be used to assess the goodness of fit of data to any distribution.

4.2 Analysis of Multiple Stocks Returns

As a majority of financial data is multivariate, it is advantageous to perform comparative analysis of multiple security returns. In some cases, one has to compare one series with another. In other cases, many variables might have to be simultaneously measured to capture the complex nature of the relationship among variables. Comparing the complexities of these factors gives the analyst a more detailed account of the relationships between selected returns, thus allowing for a better interpretation of their values and behaviors. In this section, we first compare the performance of one asset with another using graphs, then compute descriptive statistics as well as correlation matrices.

Two most commonly used graphs for comparing historical performance of more than one stock/ETF are time series plots of normalized prices and cumulative returns. The next two functions take four arguments as defined in Section  and compute normalized prices and cumulative returns.

normalizedPrices[symbols:{__String},startDate_,endDate_,period_]:=Module[
{data,mergedData,prices,dates,pricesTimeSeries,pruneprices,stock,index,result,series,newprices,indices,timeseries},
data=FinancialData[symbols,"Close",{startDate,endDate,period}];
mergedData=TimeSeriesThread[#&,data,ResamplingMethod->"NaN"];
prices=QuantityMagnitude@mergedData["Values"];
dates=mergedData["Dates"][[All,1,{1,2,3}]];
pricesTimeSeries=MapThread[Flatten[{{#1},#2},1]&,{dates,prices}];
pruneprices=Select[pricesTimeSeries[[All,{1,2}]],!MemberQ[#,"NaN"]&];
stock=pruneprices[[All,2]];
index=stock/Table[First[stock],Length[stock]];
result=MapThread[#1->{#2}&,{pruneprices[[All,1]],index}];
Do[
series=Select[pricesTimeSeries[[All,{1,i+1}]],!MemberQ[#,"NaN"]&];
newprices=series[[All,2]];
indices=newprices/Table[First[newprices],Length[newprices]];
timeseries=MapThread[#1->{#2}&,{series[[All,1]],indices}];
result=Merge[KeyUnion[{result,timeseries},Missing[]&],Identity];,
{i,2,Length[symbols]}
];
TimeSeries[Flatten[#]&/@Values[result],{Keys[result]}]
]

cumlativeReturns[symbols:{__String},startDate_,endDate_,period_]:=Module[
{data,mergedData,prices,dates,pricesTimeSeries,pruneprices,stock,index,result,series,newprices,indices,timeseries},
data=FinancialData[symbols,"Close",{startDate,endDate,period}];
mergedData=TimeSeriesThread[#&,data,ResamplingMethod->"NaN"];
prices=QuantityMagnitude@mergedData["Values"];
dates=mergedData["Dates"][[All,1,{1,2,3}]];
pricesTimeSeries=MapThread[Flatten[{{#1},#2},1]&,{dates,prices}];
pruneprices=Select[pricesTimeSeries[[All,{1,2}]],!MemberQ[#,"NaN"]&];
stock=pruneprices[[All,2]];
index=100 Accumulate[Differences[stock]/Most[stock]];
result=MapThread[#1->{#2}&,{pruneprices[[2;;,1]],index}];
Do[
series=Select[pricesTimeSeries[[All,{1,i+1}]],!MemberQ[#,"NaN"]&];
newprices=series[[All,2]];
indices=100 Accumulate[Differences[newprices]/Most[newprices]];
timeseries=MapThread[#1->{#2}&,{series[[2;;,1]],indices}];
result=Merge[KeyUnion[{result,timeseries},Missing[]&],Identity];,
{i,2,Length[symbols]}
];
TimeSeries[Flatten[#]&/@Values[result],{Keys[result]}]
]

We get normalized prices and cumulative returns and plot them for three stocks (Facebook, Inc. (FB), Costco Wholesale Corporation (COST) and Walmart Inc. (WMT)) over the period May 1, 2000, to May 30, 2019.

DateListPlot[
normalizedPrices[{"FB","COST","WMT"},{2000,5,1},{2019,5,30},"Month"],PlotLabel->"Normalized Monthly Prices Returns over Time",PlotRange->All,PlotLegends->Placed[{"FB","COST","WMT"},Below]
]

DateListPlot[cumlativeReturns[{"FB","COST","WMT"},{2010,06,10},{2019,6,30},"Month"],PlotLabel->"Cumulative Returns over Time",PlotRange->All,PlotLegends->Placed[{"FB","COST","WMT"},Below]
]

Besides graphs, we can also compute descriptive statistics and compare their performance. We define a function performanceStatisticsTable for that purpose. It takes four arguments as defined in Section  and returns a table with different types of descriptive statistics.

performanceStatisticsTable[symbols:{__String},startDate_,endDate_,period_]:=Module[
{head,n,stats,data,prices,dates,returns,var,cvar},
head={"Start Date","End Date","Holding Period Return","Arithmetic Mean","Geometric Mean","Cumulative Return","Stdandard Deviation","Minimum Return","Maximum Return","Skewness","Kurtosis","Value at Risk (VaR)","Conditional VaR (CVaR)"};
n=Length[symbols];
stats=ConstantArray[Missing[],{n,Length[head]}];
Do[
data=FinancialData[symbols[[i]],{startDate,endDate,period}];
prices=QuantityMagnitude@data["Values"];
dates=data["Dates"][[All,1,{1,2,3}]];
returns=Differences[prices]/Most[prices];
var=Sort[returns][[Floor[Length[returns] 0.05]]];
cvar=Mean@Select[returns,#≤var&];
stats[[i]]={DateString[First[dates],"Year"]<>DateString[First[dates],"Month"],DateString[Last[dates],"Year"]<>DateString[Last[dates],"Month"],100 (Last[prices]/First[prices]-1),100 Mean[returns],100 (GeometricMean[1+returns]-1),100 Total[returns],100 StandardDeviation[returns],100 Min[returns],100 Max[returns],Skewness[returns],Kurtosis[returns],100 var,100 cvar},
{i,n}
];
TableForm[Transpose[stats],TableHeadings->{head,symbols}]
]

For example, we download historical data and compute different descriptive statistics for three stocks (Walmart Inc. (WMT), Apple Inc. (AAPL) and Microsoft Corporation (MSFT) ) using monthly data over the period January 1, 2010, to March 30, 2019.

performanceStatisticsTable[{"WMT","AAPL","MSFT"},{2010,1,1},{2019,3,30},"Month"]//Text

Similarly to how we calculated an individual stock’s monthly statistics in Section ., we define the monthlyStatisticsTable function to compute monthly statistics for more than one stock given the arguments: stock ticker symbols, start date, end date and a statistical function such as Mean, StandardDeviation or Total.

statistics[symbol_String,startDate_,endDate_,function_]:=Module[
{data,prices,returns,dates,timeseries},
data=FinancialData[symbol,{startDate,endDate}];
prices=QuantityMagnitude@data["Values"];
returns=100 Differences[prices]/Most[prices];
dates=data["Dates"][[All,1,{1,2,3}]];
timeseries=Transpose[{
DateString[#,"Year"]<>DateString[#,"Month"]&/@dates[[2;;]],returns
}];
Transpose[{
DeleteDuplicates[DateString[#,"Year"]<>DateString[#,"Month"]&/@dates[[2;;]]],
(function/@GatherBy[timeseries,First])[[All,2]]
}]
]

mergeData[a_List,b_List,pad_]:=Module[
{rules,keys},
rules=Apply[#->{##2}&,{a,b},{2}];
keys=Union@@Keys@rules;
Join[List/@keys,##,2]&@@(Lookup[#,keys,pad&/@#[[1,2]]]&/@rules)
]

monthlyStatisticsTable[symbols:{__String},startDate_,endDate_,function_]:=Module[
{table,heads=Flatten[{"Date",symbols}]},
table=statistics[symbols[[1]],startDate,endDate,function];
Do[
table=mergeData[table,statistics[symbols[[i]],startDate,endDate,function],Missing[]];,{i,2,Length[symbols]}
];
Dataset[AssociationThread[heads->#]&/@table]
]

For example, we compute the monthly cumulative returns for four stocks (Walmart Inc. (WMT), Apple Inc. (AAPL), Microsoft Corporation (MSFT) and Netflix, Inc. (NFLX)) over the period January 1, 2010, to June 30, 2019, and create a dataset. The first column represents year and month, the first four digits for the year and the last two digits for the month.

monthlyStatisticsTable[{"WMT","AAPL","MSFT","NFLX"},{2010,01,1},{2019,06,30},Total]

Similarly, box-and-whisker charts, paired histograms, paired smooth histograms and matrix scatterplots are often used to examine multivariate data. The BoxWhiskerChart function can be used to make a box plot that gives a glimpse of the distribution of the given dataset. You can see the statistical information by hovering over the boxes in the plot. The PairedSmoothHistogram and SmoothKernelDistribution functions are used to create paired histogram and smooth distribution plots. They can be used to compare how two datasets are distributed. The PairwiseScatterPlot function from the Statistical Plots package can be used to make scatter plots of multivariate data. It creates scatter plots comparing the data in each column against other columns. More complex analysis of multivariate data can be done using functions from the Multivariate Statistics package. The package contains functions to compute descriptive statistics for multivariate data and distributions derived from the multivariate normal distribution. All these functions are well explained in the official documentation.

5. Rolling-Window Performance Analysis

Rolling-window performance analysis is a simple technique to access variability of the statistical performance measures. For example, if we want to access the stability of mean or standard deviation of returns on a stock over time, we can choose rolling window (the number of consecutive observations per rolling window), estimate the mean or standard deviation and plot series of the estimates. A little fluctuation is normal but large fluctuations indicate a shift in the values of the estimate. Built-in functions such as MovingAverage, MovingMap and TimeSeriesAggregate are useful for rolling-window performance analysis. In this section, we are going to show a few examples of how to compute rolling-window-based performance statistics.

We define plotRollingWindowStatistics that can be used to plot rolling-window statistics given its five inputs: stock ticker symbol, start date, end date, size of window in days and function to apply, which can be any built-in or user-defined function.

plotRollingWindowStatistics[symbol_String,startDate_,endDate_,windowInDays_,function_]:=
DateListPlot[
MovingMap[
function,
FinancialData[symbol,"FractionalChange",{startDate,endDate}],
Quantity[windowInDays,"Days"]
],
PlotTheme->"Detailed"
]

For example, we compute and plot the 90-day rolling mean to standard deviation ratio on Walmart's daily stock returns over the period January 1, 2001, to March 30, 2019.

plotRollingWindowStatistics["WMT",{2001,01,01},{2019,03,30},90,Mean[#]/StandardDeviation[#]&]

Similarly, we define the function plotRollingCorrelation to compute the rolling correlation of two series and apply it together with MovingMap to the desired data.

plotRollingCorrelation[{symbol1_String,symbol2_String},startDate_,endDate_,windowInDays_]:=Module[
{returns,correlationFunction},
returns=TimeSeriesThread[#&,FinancialData[{symbol1,symbol2},"FractionalChange",{startDate,endDate}],ResamplingMethod->Missing[]];
correlationFunction=N[Apply[Correlation,Transpose[#]]]&;
DateListPlot[MovingMap[correlationFunction,returns,Quantity[windowInDays,"Days"]],PlotRange->All,PlotTheme->"Detailed"]
]

We use the plotRollingCorrelation function to plot the 90-day rolling correlation of daily returns on two stocks, WMT and COST, for the period from March 30, 2009, to March 30, 2019.

plotRollingCorrelation[{"WMT","COST"},{2009,01,30},{2019,03,30},90]

Sometimes, it is also useful to store these time-varying descriptive statistics as a dataset so that we can use them in the subsequent analysis. The function rollingStatisticsTable, given its input, computes the geometric mean, standard deviation and the ratio of the arithmetic mean to the standard deviation on a rolling-window basis.

rollingStatisticsTable[symbol_String,startDate_,endDate_,windowInDays_]:=Module[
{statisticsToCompute,statisticsComputed,statisticsValues,dates},
statisticsToCompute=Flatten[{100 (GeometricMean[1+#]-1),StandardDeviation[#],Mean[#]/StandardDeviation[#]}]&;
statisticsComputed=MovingMap[
statisticsToCompute,
FinancialData[symbol,"FractionalChange",{startDate,endDate}],
Quantity[windowInDays,"Days"]
];
statisticsValues=QuantityMagnitude[statisticsComputed]["Values"];
dates=DateObject[#]&/@statisticsComputed["Dates"][[All,1,{1,2,3}]];
Dataset[AssociationThread[{"Date","GM","Std. Dev.","AM / Std. Dev."}->#]&/@Transpose@Insert[Transpose[statisticsValues],dates,1]]
]

For example, we compute the 90-day rolling-window geometric mean (GM), standard deviation (Std. Dev.) and arithmetic mean to standard deviation ratio (AM/Std. Dev.) of Walmart’s daily stock returns over the period July 1, 2018, to October 30, 2019. You can scroll through the dataset.

rollingStatisticsTable["WMT",{2018,07,01},{2019,10,30},90]

6. Portfolio Optimization

Mean-variance analysis is one of the foundations of financial economics. Portfolio optimization is essential, whether it be in professional or personal financial planning. In this section, we are going to show how to implement the most commonly used optimization techniques in finance using historical returns. We want to point out that future returns on investment depend on expected returns and other conditioning information, not on the past returns. Past returns are used only for illustration and do not guarantee future returns.

Define the following variables:

• rf is the risk-free rate

• (Subscript[w, i]) is the proportion of wealth invested in security i

• (Subscript[μ, i]) is the average return on security i

• (Subscript[σ, i]^2) is the variance of security i

• (Subscript[σ, i,j]) is the covariance between securities i and j

• (Subscript[ρ, i,j]) is the correlation between securities i and j

Then we define the vectors of mean returns and weights and the covariance matrix:

μ=(({{(Subscript[μ, 1])}, {⋮}, {(Subscript[μ, n])}})), w=(({{(Subscript[w, 1])}, {⋮}, {(Subscript[w, 2])}})) and Σ=(({{(Subscript[σ, 11]), (Subscript[σ, 12]), …, (Subscript[σ, 1,n])}, {(Subscript[σ, 21]), (Subscript[σ, 22]), …, (Subscript[σ, 2,n])}, {⋮, ⋮, ⋱, ⋮}, {(Subscript[σ, n,1]), (Subscript[σ, n,2]), …, (Subscript[σ, n,n])}})).

The formulas for the portfolio mean and variance are (Subscript[μ, p])=μ·w=(∑_( i=1)^n)(Subscript[μ, i])(Subscript[w, i]) and (Subscript[σ, p]^2)=w·Σ·w=(∑_( i=1)^n)(∑_( j=1)^n)(Subscript[w, i])(Subscript[w, j])(Subscript[σ, i,j]), respectively. The corresponding Mathematica code is μ.w and w.Σ.w.

In order to compute portfolio statistics, we need returns data. We can use the getReturnsData function to download historical returns data. It takes four arguments as defined in Section  and gives a matrix of returns. Most functions in this section use the getReturnsData function, so run it before you run other functions.

getReturnsData[symbols:{__String},startDate_,endDate_,period_]:=Module[
{data,mergedData,prices,pruneprices},
data=FinancialData[#,{startDate,endDate,period}]&/@symbols;mergedData=TimeSeriesThread[#&,data,ResamplingMethod->Missing[]];prices=QuantityMagnitude[mergedData["Values"]];
pruneprices=Select[prices,! MemberQ[#,QuantityMagnitude[Missing[]]]&];
Differences[pruneprices]/Most[pruneprices]
]

To compute basic portfolio statistics such as portfolio mean, variance, standard deviation and Sharpe ratio, we can use portfolioStatistics, which takes six arguments. The first four arguments are as defined in Section  and the other two are a list of weights and the optional risk-free rate.

portfolioStatistics[symbols:{__String},startDate_,endDate_,period_,weights_List,riskfreeRate_:0]:=Module[
{returns=getReturnsData[symbols,startDate,endDate,period],portfolioMean,portfolioVariance},
portfolioMean=weights.Mean@returns;
portfolioVariance=weights.Covariance@returns.weights;
TableForm[{
{
"Portfolio\nMean",
"Portfolio\nVariance",
"Portfolio\nStd. Dev.",
"Sharpe\nRatio"
},
{
portfolioMean,
portfolioVariance,
Sqrt[portfolioVariance],
(portfolioMean-riskfreeRate)/Sqrt[portfolioVariance]
}
},TableAlignments->Center]
]

For example, we compute the portfolio mean, variance, standard deviation and Sharpe ratio for the portfolio that consists of the stock returns of five companies: Apple (AAPL), Walmart (WMT), Boeing (BA), 3M and Exxon Mobil (XOM), using monthly returns over the period January 1, 2009, to May 30, 2019.

portfolioStatistics[{"AAPL","WMT","BA","MMM","XOM"},{2009,1,1},{2019,5,30},"Month",{0.2,0.2,0.2,0.2,0.2},0.0023]//Text

The plotFrontier function plots the Markowitz portfolio frontier; it takes a matrix of returns obtained from getReturnsData as its only argument. The function uses the concept that any two efficient portfolios are enough to establish the whole portfolio frontier, as first proved by Black []. It accepts any ListPlot option.

plotFrontier[returns_List,options:OptionsPattern[ListPlot]]:=Module[
{μ,σ,ones,wt1,wt2,w1,means,standardDeviations},
μ=Mean[returns];
σ=Covariance[returns];
ones=ConstantArray[1,Length[μ]];
wt1=Inverse[σ].ones/ones.Inverse[σ].ones;
wt2=Inverse[σ].μ/ones.Inverse[σ].μ;
means=w1 wt1.μ+(1-w1) wt2.μ/.w1->Range[-0.75,1.75,0.01];
standardDeviations=Sqrt[w1^2 wt1.σ.wt1+(1-w1)^2 wt2.σ.wt2+2 w1 (1-w1) wt2.σ.wt1]/.w1->Range[-0.75,1.75,0.01];
Labeled[
ListLinePlot[Transpose[{standardDeviations,means}],PlotTheme->"Detailed",AxesOrigin->{Min[standardDeviations] 0.75,Min[means] 1.5},options],Text/@{"Portfolio Standard Deviation((Subscript[σ, p]))","Portfolio Mean ((Subscript[μ, p]))"},
{Bottom,Left},RotateLabel->True
]
]

For example, we plot the portfolio frontier for the portfolio that consists of stock returns of five companies: Apple (AAPL), Walmart (WMT), Boeing (BA), 3M and Exxon Mobil (XOM), using monthly returns over the period January 1, 2009, to May 30, 2019.

plotFrontier[
getReturnsData[{"AAPL","WMT","BA","MMM","XOM"},{2009,1,1},{2019,5,30},"Month"],
PlotLabel->"Portfolio Frontier",ColorFunction->"DarkRainbow"]

Next, we solve two kinds of portfolio problems: global minimum variance portfolio and tangency portfolio. In terms of the notation defined earlier in this section, the global minimum variance portfolio can be obtained by minimizing (Subscript[σ, p]^2) subject to (∑_(i=1)^n)(Subscript[w, i])=1 and solving for (Subscript[w, i]). Its solution can be obtained with the built-in Mathematica function NMinimize.

The function globalMininimumVariancePortfolioWeights computes weights, returning the portfolio allocation on stocks considered for a global minimum variance portfolio.

globalMininimumVariancePortfolioWeights[symbols:{__String},startDate_,endDate_,period_]:=Module[
{returns,weights},
returns=getReturnsData[symbols,startDate,endDate,period];
weights=Subscript[w,#]&/@Range[Dimensions[returns][[2]]];
weights /. NMinimize[
{
weights.(12Covariance@returns).weights,
Total[weights]==1
},
weights][[2]
]
]

We compute the global minimum variance portfolio weights using the monthly stock returns of five companies: Apple (AAPL), Walmart (WMT), Boeing (BA), 3M and Exxon Mobil (XOM), over the period January 1, 2009, to May 30, 2019.

globalMininimumVariancePortfolioWeights[{"AAPL","WMT","BA","MMM","XOM"},{2009,1,1},{2019,5,30},"Month"]

Similarly, the tangency portfolio can be obtained by maximizing ((Subscript[μ, p])-rf)/(Subscript[σ, p]), where rf is the risk-free rate (a constant in this case), subject to (∑_(i=1)^n)(Subscript[w, i])=1 and solving for (Subscript[w, i]). The solution uses the built-in Mathematica function NMaximize. The tangencyPortfolioWeights function computes the tangency portfolio weights given its five inputs, four as defined in the Section  and the risk-free rate.

tangencyPortfolioWeights[symbols:{__String},startDate_,endDate_,period_,rf_]:=Module[
{returns,weights},
returns=getReturnsData[symbols,startDate,endDate,period];
weights=Subscript[w,#]&/@Range[Dimensions[returns][[2]]];
weights /. NMaximize[
{
(weights.(12 Mean@returns)-rf)/Sqrt[weights.(12 Covariance@returns).weights],
Total[weights]==1
},
weights][[2]]
]

Assuming a monthly risk-free rate of 0.1667 percent and using monthly data over the period January 1, 2009, to May 30, 2019, we calculate the tangency portfolio weights for our portfolio of five stocks: (Apple (AAPL), Walmart (WMT), Boeing (BA), 3M and Exxon Mobil (XOM)).

tangencyPortfolioWeights[{"AAPL","WMT","BA","MMM","XOM"},{2009,1,1},{2019,5,30},"Month",0.001667]

Portfolio optimization using the Wolfram Language is very flexible. We can formulate any kind of portfolio and use built-in functions such as NMinimize, NMaximize or QuadraticOptimization to get numerical solutions to the portfolio problem.

7. Portfolio Risk Decomposition

In this section, we concentrate on how to decompose a measure of portfolio risk (portfolio standard deviation) into risk contribution from individual assets included in the portfolio. It helps to see how individual assets influence portfolio risk. When risk is measured by standard deviation, we can use Euler’s theorem for decomposing risk into asset-specific risk contribution. Euler’s theorem provides an additive decomposition of a homogeneous function. For reference, see Campolieti and Makarov []. Using Euler’s theorem, we can define the percentage contribution to portfolio standard deviation of an asset as mc×wt/(Subscript[σ, p]), where mc is the marginal contribution of the asset, wt is the weight of the asset, and (Subscript[σ, p]) is the portfolio standard deviation.

We define the function decomposePortRisk for portfolio risk decomposition. It takes five arguments, four arguments as defined in Section  and a list of portfolio weights; it returns a bar chart representing the individual asset's risk contribution to the portfolio standard deviation.

decomposePortRisk[symbols:{__String},startDate_,endDate_,period_,wt_List]:=Module[
{data,mergedData,prices,pruneprices,returns,Σ,σ,mc},
data=FinancialData[#,{startDate,endDate,period}]&/@symbols;
mergedData=TimeSeriesThread[#&,data,ResamplingMethod->Missing[]];
prices=QuantityMagnitude[mergedData["Values"]];
pruneprices=Select[prices,!MemberQ[#,QuantityMagnitude[Missing[]]]&];
returns=Differences[pruneprices]/Most[pruneprices];
Σ=12Covariance@returns;
Labeled[BarChart[wt (wt.Σ)/(wt.Σ.wt),ChartLabels->symbols,ChartStyle->"DarkRainbow"],Text/@{"Asset Symbols","Percent Contribution to Risk"},{Bottom,Left},RotateLabel->True]
]

We calculate the risk contribution of each asset in a portfolio that consists of five stocks using the historical monthly returns over the period January 30, 2010, to May 30, 2019 and make a bar chart.

decomposePortRisk[{"AAPL","WMT","BA","MMM","XOM"},{2010,01,30},{2019,05,30},"Month",{0.20,0.20,0.2,0.2,0.2}]

8. Factor Models

Currently, factor models are widely accepted and used in finance to construct portfolios, to evaluate portfolio performance and for risk analysis. Factor models are regression models. We can use the built-in function LinearModelFit to estimate and evaluate the appropriateness of the regression models. In addition, we can download all factor data directly from Prof. Kenneth French’s data library. Before we apply factor models to real-world data, we need data in the form

(I,D),

where

I=(({{(Subscript[f, 11]), (Subscript[f, 12]), …, (Subscript[f, 1,k])}, {(Subscript[f, 21]), (Subscript[f, 22]), …, (Subscript[f, 2,k])}, {⋮, ⋮, ⋱, ⋮}, {(Subscript[f, n,1]), (Subscript[f, n,2]), …, (Subscript[f, n,k])}}))

is a matrix of values of independent variables and

D=(({{(Subscript[y, 1])}, {(Subscript[y, 2])}, {⋮}, {(Subscript[y, n])}}))

is a vector of values of a dependent variable.

Commonly used factor models are summarized in Table . We can find more about the factor models in Fama and French []. We use the following notation: r-rf for excess return on the security or portfolio, MKT for excess return on the value-weighted market portfolio, SML for return on a diversified portfolio of small-capitalization stocks minus the return on a diversified portfolio of large-capitalization stocks, HML for difference in the returns on diversified portfolios of high-book-to-market stocks and low-book-to-market stocks, MOM for difference in returns on diversified portfolios of the prior year’s winners and losers, RMW for difference between the returns on diversified portfolios of stocks with robust and weak profitability, CMA for difference between the returns on diversified portfolios of the stocks of low and high investment firms, α for risk-adjusted return on security or portfolio and (Subscript[β, MKT]), (Subscript[β, SMB]), (Subscript[β, HML]), (Subscript[β, MOM]), (Subscript[β, RMW]) and (Subscript[β, CMA]) for betas or factor loadings.

Before we estimate these models using real data, we need factors data. We download five factors and the momentum factor data from Kenneth French’s website (http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html) and define variables to store them.

FamaFrench5FactorData=First@Import["http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Research_Data_5_Factors_2x3_CSV.zip","*"];

FamaFrenchMomentumFactor=First[Import["http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Momentum_Factor_CSV.zip","*"]];

The function factorsData takes only two arguments, the start date and end date, and returns time series of all factors data. The start date and end date must be specified as date objects.

factorsData[startDate_,endDate_]:=Module[
{startDateFamaFrenchFactor,startDateMomentumFactor,lastDate,dateFactorAssociation,keys,mergedFactorData,dateRange,mergedFactorDataFinal},
startDateFamaFrenchFactor={1963,07,1};
startDateMomentumFactor={1927,01,1};
lastDate=DatePlus[Today,{-4,"Month"}];
dateFactorAssociation=Apply[#->{##2}&,{FamaFrench5FactorData[[5;;Length[DateRange[startDateFamaFrenchFactor,lastDate,"Month"][[All,1;;2]]]+5]],FamaFrenchMomentumFactor[[15;;Length[DateRange[startDateMomentumFactor,lastDate,"Month"][[All,1;;2]]]+15]]},{2}];
keys=Union@@Keys@dateFactorAssociation;
mergedFactorData=Join[List/@keys,##,2]&@@(Lookup[#,keys,Missing[]&/@#[[1,2]]]&/@dateFactorAssociation);
dateRange=DateRange[startDate,endDate,"Month"][[All,1;;3]];
mergedFactorDataFinal=Select[mergedFactorData,(startDate[[1]] 100+startDate[[2]])≤#[[1]]≤(endDate[[1]] 100+endDate[[2]])&];
Prepend[
Table[
Flatten[{{dateRange[[i,1;;2]]},mergedFactorDataFinal[[i,{2,3,4,5,6,8,7}]]},1],
{i,Length[mergedFactorDataFinal]}
],
{"Date","Mkt-RF","SMB","HML","RMW","CMA","Mom","RF"}
]
]

Similarly, the monthlyStockReturns function can be used to get a stock’s monthly returns data. It takes three arguments, a ticker symbol of any publicly traded company, and the start and end dates for the analysis period.

monthlyStockReturns[symbol_String,startDate_,endDate_]:=Module[
{data,prices,returns,dates},
data=FinancialData[symbol,"Close",{startDate,endDate,"Month"}];
prices=QuantityMagnitude[data["Values"]];
returns=100 Differences[prices]/Most[prices];
dates=data["Dates"][[2;;,1,{1,2}]];
Prepend[Transpose[{dates,returns}],{"Date",symbol}]
]

Using factorsData and monthlyStockReturns, we define dataforNFactorModel to combine factors and returns data. The dataforNFactorModel takes four arguments (the symbol of the stock for which we want to estimate the factor model, the start date, the end date and an integer that represents the number of factors) and returns a data matrix suitable for LinearModelFit.

dataforNFactorModel[symbol_String,startDate_,endDate_,nfac_/;MemberQ[{1,3,4,5},nfac]]:=Module[
{returns,dateAssociationData,mergedData,mergedDatafinal,columnsToInclude},
returns=monthlyStockReturns[symbol,startDate,endDate];
dateAssociationData=Apply[#->{##2}&,{factorsData[startDate,endDate][[2;;]],returns[[2;;]]},{2}];
mergedData=Merge[KeyUnion[dateAssociationData,Missing[]&],Identity];
mergedDatafinal=Select[Flatten[#]&/@Values[mergedData],!MemberQ[#,Missing[]]&];
columnsToInclude=Switch[nfac,1,{1},3,{1,2,3},4,{1,2,3,6},5,{1,2,3,4,5}];
Table[
Flatten[{mergedDatafinal[[i,columnsToInclude]],mergedDatafinal[[i,8]]-mergedDatafinal[[i,7]]}],
{i,Length[mergedDatafinal]}
]
]

Next, we estimate different factor models for Apple’s stock using monthly data from October 1, 2008, to March 30, 2019. We estimate the capital asset pricing model (CAPM) with market factor (MKT).

capm=LinearModelFit[dataforNFactorModel["AAPL",{2008,10,1},{2019,3,30},1],{MKT},{MKT}];
capm["ParameterTable"]

Grid[Transpose[{#,capm[#]}&[{"AdjustedRSquared","AIC","BIC","RSquared","DurbinWatsonD"}]],Alignment->Left]//Text

We estimate the Fama–French three-factor model with market, size and value factors (MKT, SMB, HML).

threeFactorModel=LinearModelFit[dataforNFactorModel["AAPL",{2008,10,1},{2019,3,30},3],{MKT,SMB,HML},{MKT,SMB,HML}];
threeFactorModel["ParameterTable"]

Grid[Transpose[{#,threeFactorModel[#]}&[{"AdjustedRSquared","AIC","BIC","RSquared","DurbinWatsonD"}]],Alignment->Left]//Text

Similarly, we estimate the Carhart four-factor model with market, size, value and momentum factors (MKT, SMB, HML, MOM).

fourFactorModel=LinearModelFit[dataforNFactorModel["AAPL",{2008,10,1},{2019,3,30},4],{MKT,SMB,HML,MOM},{MKT,SMB,HML,MOM}];
fourFactorModel["ParameterTable"]

Grid[Transpose[{#,fourFactorModel[#]}&[{"AdjustedRSquared","AIC","BIC","RSquared","DurbinWatsonD"}]],Alignment->Left]//Text

Finally, the Fama–French five-factor model with market, size, value, profitability and investment factors (MKT, SMB, HML, RMW, CMA) is estimated as follows.

fiveFactorModel=LinearModelFit[dataforNFactorModel["AAPL",{2008,10,1},{2019,3,30},5],{MKT,SMB,HML,RMW,CMA},{MKT,SMB,HML,RMW,CMA}];
fiveFactorModel["ParameterTable"]

Grid[Transpose[{#,fiveFactorModel[#]}&[{"AdjustedRSquared","AIC","BIC","RSquared","DurbinWatsonD"}]],Alignment->Left]//Text

Once the model is estimated, we must access different properties related to data and fitted models. To assess how well the model fits the data and how well the model meets the assumptions, there are many built-in functions.To learn more about obtaining diagnostic information, see the properties of LinearModelFit.

9. Stock Portfolio Performance Measures

There are various measures to evaluate the performance and risk of portfolios. Most of these measures are used to evaluate a portfolio of interest against a chosen benchmark, by taking a snapshot of the past or considering the entire historical picture. We will compute some common metrics often employed by investors while analyzing performance measures. All the computations are based on the formulas developed in Bacon []. Most common measures are summarized in the function computePortfolioPerformanceTable. The function takes four arguments:

• ticker symbols for the test assets and the benchmark asset

• periodic risk-free rate

• time period

• frequency of data

We repeat the definition of the getReturnsData function defined earlier to make this section self-contained.

getReturnsData[symbols_:{__String},startDate_,endDate_,period_]:=Module[
{data,mergedData,prices,pruneprices},
data=FinancialData[#,{startDate,endDate,period}]&/@symbols;
mergedData=TimeSeriesThread[#&,data,ResamplingMethod->Missing[]];prices=QuantityMagnitude[mergedData["Values"]];
pruneprices=Select[prices,! MemberQ[#,QuantityMagnitude[Missing[]]]&];
Differences[pruneprices]/Most[pruneprices]
]

computePortfolioPerformanceTable[testAssets_:{__String},benchmarkAsset_String,riskfreeRate_,startDate_,endDate_,period_]:=Module[{stockSymbols,returns,benchmarkReturns,totalReturns,n,mean,standardDeviation,skewness,kurtosis,SharpeRatio,adjustedSharpeRatio,trackingError,informationRatio,minReturns,SortinoRatio,historicalVAR,historicalCVAR,quantile,normalVAR,normalCVAR,mSquare,performanceMeasures,head},
stockSymbols=Flatten[{benchmarkAsset,testAssets}];
returns=getReturnsData[stockSymbols,startDate,endDate,period];
benchmarkReturns=returns[[All,1]];
performanceMeasures=ConstantArray["NaN",{Length[stockSymbols[[2;;]]],19}];
Do[
totalReturns=returns[[All,i]];
n=Length[benchmarkReturns];
mean=Mean[totalReturns];
standardDeviation=StandardDeviation[totalReturns];
skewness=n/((n-1) (n-2)) Total[(totalReturns-mean)^3]/standardDeviation^3;
kurtosis=n (n+1)/((n-1) (n-2) (n-3)) Total[(totalReturns-mean)^4]/standardDeviation^4-(3 (n-1)^2)/((n-2) (n-3));
SharpeRatio=(mean-riskfreeRate)/standardDeviation;
adjustedSharpeRatio=SharpeRatio (1+skewness/6 SharpeRatio-kurtosis/24 SharpeRatio^2);
trackingError=StandardDeviation[totalReturns-benchmarkReturns];
informationRatio=(mean-Mean[benchmarkReturns])/trackingError;
minReturns=Min[#-riskfreeRate,0]&/@totalReturns;
SortinoRatio=(mean-riskfreeRate)/Sqrt[Total[minReturns^2]/(Length[minReturns]-1)];
historicalVAR=Sort[totalReturns][[Floor[Length[totalReturns] 0.05]]];
historicalCVAR=Mean[Select[totalReturns,#≤historicalVAR&]];
quantile=Quantile[NormalDistribution[],0.05];
normalVAR=mean+standardDeviation quantile;
normalCVAR=mean-standardDeviation (PDF[NormalDistribution[],quantile]/CDF[NormalDistribution[],quantile]);
mSquare=((mean-riskfreeRate)/standardDeviation-(Mean[benchmarkReturns]-riskfreeRate)/StandardDeviation[benchmarkReturns]) StandardDeviation[benchmarkReturns];
performanceMeasures[[i-1]]={
mean,
standardDeviation,
100 (GeometricMean[1+totalReturns/100]-1),
Total[totalReturns],
skewness,
kurtosis,
SharpeRatio,
adjustedSharpeRatio,
trackingError,
informationRatio,informationRatio (1+skewness/6 informationRatio-kurtosis/24 informationRatio^2),
SortinoRatio,
Min[Accumulate/@SplitBy[totalReturns,Sign]],
Sqrt[Total[(Boole@Thread[totalReturns<0] totalReturns)^2]/n],
historicalVAR,
normalVAR,
historicalCVAR,
normalCVAR,
mSquare
},
{i,2,Length[stockSymbols]}
];
head={
"Arithmetic Mean",
"Standard Deviation",
"Geometric Mean",
"Cumulative Returns",
"Sample Skewness",
"Sample Excess Kurtosis",
"Sharpe Ratio",
"Adjusted Sharpe Ratio",
"Tracking Error",
"Information Ratio",
"Adjusted Information Ratio",
"Sortino Ratio",
"Continuous Drawdown",
"Pure Downside Risk",
"Historical 95% VaR",
"Normal 95% VaR",
"Historical 95% CVaR",
"Normal 95% CVaR",
"M-Square"
};
Prepend[
Transpose@Prepend[performanceMeasures,head],Join[{"Measures"},stockSymbols[[2;;]]]
]//TableForm]

The next example uses stocks, although these measures are used to evaluate the performance of portfolios, mutual funds and exchange-traded funds. We evaluate the performance of (Walmart Inc. (WMT), Apple Inc. (AAPL) and Microsoft Corporation (MSFT) ) against the S&P 500 index (^SPX) using 0.0016 as a monthly risk-free rate and month as the data frequency over the period January 1, 1995, to March 30, 2019.

computePortfolioPerformanceTable[{"WMT","AAPL","MSFT"},"^SPX",0.0016,
{1995,01,01},{2019,03,30},"Month"]//Text

10. Interactive Graphics and Technical Analysis of Stock Prices

Short-term traders commonly use interactive graphics and technical indicators of stock prices to profit from stocks that may be overbought or oversold. Much is based on market sentiment, but also market timing. When a stock is oversold, the price is low and people want to buy. In comparison, when a stock is overbought, the price is between its normal range and higher, and people do not want to buy, or may want to short sell. Many technical indicators are used to determine a given stock’s peak or bottom price and how to take advantage of that information. Three of the most useful functions for technical analysis of stocks are TradingChart, InteractiveTradingChart and CandlestickChart. The documentation provides a comprehensive set of examples on how to use them.

We show one example of how to use the TradingChart function and one example of how to use InteractiveTradingChart.

You can choose the chart type and from over 100 technical indicators, which are divided into eight groups:

• basic

• moving average

• market strength

• momentum

• trend

• volatility

• band

• statistical

FinancialIndicator[All] displays all the available technical indicators.

FinancialIndicator[All]

Here is the basic format:

TradingChart[{symbol,{startDate, endDate},{indicator1,indicator2,…}]

Alternatively, use:

TradingChart[timeseriesdata,{indicator1,indicator2,…}]

The time series data must be of the form:

{
{date1,{open1,high1,low1,close1,volume1}},
{date2,{open2,high2,low2,close2,volume2}},
…
}

The historical open, high, low, close and volume data retrieved from FinancialData can also be used as data input. The TradingChart function has many options that can be used to enhance graphics. We produce a chart using historical prices of Apple's stock and volume over the period January 1, 2018, to March 30, 2019.

The top of the chart shows a plot of historical prices and 50-day and 200-day moving averages. The second part shows the historical volume. The last two parts show plots of two indicators, the commodity channel index (CommodityChannelIndex) and the relative strength index (RelativeStrengthIndex).

A good introduction to technical indicators can be found in standard references, including at Fidelity Learning Center [].

TradingChart[{"AAPL",{{2018,01,01},{2019,03,30}}},{"Volume",FinancialIndicator["SimpleMovingAverage",50],FinancialIndicator["SimpleMovingAverage",200],"CommodityChannelIndex","RelativeStrengthIndex"}]

The InteractiveTradingChart function provides a point-and-click interactive chart, with a similar setup:

InteractiveTradingChart[{symbol,{startDate, endDate},{indicator1,indicator2,…}]

Alternatively, use:

InteractiveTradingChart[timeseriesdata,{indicator1,indicator2,…}]

For example, we make a chart showing prices, volume and indicators for historical data of Apple's stock over the period January 1, 2018, to March 30, 2019. The InteractiveTradingChart function provides a user-friendly environment where you can drag a slider to view different parts of the chart or you can choose different indicators with point-and-click.

InteractiveTradingChart[{"AAPL",{{2018,01,01},{2019,03,30}}},{"Volume",FinancialIndicator["SimpleMovingAverage",50],FinancialIndicator["SimpleMovingAverage",200],"CommodityChannelIndex","RelativeStrengthIndex"}]

11. Essentials of Bond Mathematics

A bond is a long-term debt instrument in which a borrower agrees to make payments of principal and interest, on specific dates, to the holders of the bond. When it comes to analysis and pricing of bonds and computing returns, convexity and duration are important concepts. When a bond is traded between coupon payment dates, its price has two parts: the quoted price and the accrued interest. The quoted price is net of accrued interest, and does not include accrued interest. Accrued interest is the proportion or share of the next coupon payment. The full price is the price of a bond including any interest that has accrued since issue of the most recent coupon payment. Similarly, yield to maturity is the rate of return earned on a bond if it is held to maturity. Duration is a measure of the average length of time for which money is invested in a coupon bond. Convexity estimates the change in the bond price given a change in the yield to maturity, assuming a nonlinear relationship between the two. The built-in functions FinancialBond and FindRoot can be used to compute various properties including value of the bond, accrued interest, yield, duration, modified duration, convexity, and so on. This section provides a few examples of how to use FinancialBond and how the concepts of bond convexity and duration can be used in bond portfolio management.

We discuss zero-coupon bonds first. The zero-coupon bond does not make coupon payments. The only cash payment is the face value of the bond on the maturity date. The yield to maturity (ytm) for a zero-coupon bond with n periods to maturity, current price P and face value FV can be obtained by solving P=FV/(((1+ytm))^n).

For example, we compute the yield to maturity of a zero-coupon bond with a $10,000 face value, time to maturity 4 years and current price $9,662 using NSolve and FindRoot.

NSolve[9662==(10000/(((1+ytm))^4))&&ytm≥0,ytm,Reals]

FindRoot[
FinancialBond[
{"FaceValue"->10000,"Coupon"->0,"Maturity"->4,"CouponInterval"->1},
{"InterestRate"->ytm,"Settlement"->0}
]==9662,
{ytm,.1}
]

Similarly, FindRoot can also be used to compute the yield to maturity of a nonzero coupon payment bond.

For example, we compute the yield to maturity of a $1,000 par value 10-year bond with 5% semiannual coupons issued on June 20, 2013, with maturity date of June 20, 2023, selling for $920 on September 15, 2018.

FindRoot[
FinancialBond[
{"FaceValue"->1000,"Coupon"->0.05,"Maturity"->{2023,6,20},"CouponInterval"->(1/2)},
{"InterestRate"->y,"Settlement"->{2018,9,15}}
]==920,
{y,.1}
]

FinancialBond can also be used to compute the price, duration, modified duration and convexity of a bond. For example, we compute those values for a bond with 8% yield, 8% annual coupons, 10-year maturity and $1,000 face value.

FinancialBond[
{"FaceValue"->1000,"Coupon"->0.08,"Maturity"->10,"CouponInterval"->1},
{"InterestRate"->.08,"Settlement"->0},
{"Value","Duration","ModifiedDuration","Convexity"}
]

There are different approaches to bond portfolio management. We concentrate here on a liability-driven portfolio strategy, in which the characteristics of the bonds that are held in the portfolio are coordinated with those of the liabilities the investor is obligated to pay. The matching techniques can range from an attempt to exactly match the levels and timing of the required cash payments to more general approaches that focus on other investment characteristics, such as setting the average duration or convexity of the bond portfolio equal to that of the underlying liabilities. One specific example would be to construct the portfolio so that the duration of the bond portfolio is equal to the duration of cash obligation and the total money invested in the bond portfolio today is equal to the present value of the future cash obligations.

To illustrate the concept of bond portfolio management, assume that we have an obligation to pay $1,000,000 in 10 years and there are two bonds available for investment. The first bond matures in 30 years with $100 face value and annual coupon payment of 6%. The second bond matures in 10 years with $100 face value and annual coupon payment of 5%. The yield to maturity is 9% on both bonds. We can decide on how much to invest in each bond so that the overall portfolio is immunized against changes in the interest rate.

We compute the duration of each bond using FinancialBond, which gives that the duration of bond 1 is 11.88 and that of bond 2 is 6.75. Assuming that the proportion of money invested in bonds 1 and 2 is (Subscript[w, 1]) and (Subscript[w, 2]), the immunized portfolio is found by solving the simultaneous equations:

{({{(Subscript[w, 1])+(Subscript[w, 2])=1}, {11.88(Subscript[w, 1])+6.75(Subscript[w, 2])=10}})

These two equations can be solved using Solve.

Solve[(Subscript[w, 1])+(Subscript[w, 2])==1&&11.88 (Subscript[w, 1])+6.75 (Subscript[w, 2])==10,{(Subscript[w, 1]),(Subscript[w, 2])}]

The result shows how much money should be allocated to each bond. More examples can be found in Benninga []. A more general approach to bond portfolio management can be solved by using linear programming. It is beyond the scope of this article to introduce linear programming.

12. Binomial and Black–Scholes–Morton Stock Option Pricing Models

The most popular options pricing models are the binomial model and the Black–Scholes–Morton option pricing formulas for European options. In the next two subsections, we discuss these models and their implementation.

12.1 Binomial Option Pricing Model

Following the notation from Hull [], define:

(Subscript[S, 0]): the current stock price

K: the strike price

rf: the annual risk-free rate

σ: the annual standard deviation of the stock returns

T: the time to expiration of the option in years

n: the total number of up and down moves

j: the number of upward moves.

Then n-j is the number of down moves on the tree.

In terms of these variables, we can define:

• time per period Δ t=T/n

• up factor u=(e^(σ(Sqrt[Δ t])))

• down factor d=(e^(-σ(Sqrt[Δ t])))

• probability of an up move p=(((e^(rf Δ t))-d)/((u-d) (e^(rf Δ t))))

• probability of a down move q=1-p

• stock price (Subscript[S, 0])(u^i)(d^(n-j)) at node n

• payoff from a European call max[(Subscript[S, 0])(u^j)(d^(n-j))-K,0]

• payoff from a European put max[K-(Subscript[S, 0])(u^j)(d^(n-j)),0]

In the risk-neutral world, the price of the call and put using the n-period binomial options pricing model can be computed as:

call=(e^(-rf T))(∑_(j=0)^n)(({{n}, {j}}))(p^j)(q^(n-j))max[(Subscript[S, 0])(u^j)(d^(n-j))-K,0],

put=(e^(-rf T))(∑_(j=0)^n)(({{n}, {j}}))(p^j)(q^(n-j))max[K-(Subscript[S, 0])(u^j)(d^(n-j)),0].

The functions binomialCallPrice and binomialPutPrice calculate the prices of European call and put options; they output the option price. Each function takes six arguments as defined at the beginning of this subsection.

binomialCallPrice[S0_,K_,σ_,rf_,T_,n_]:=Module[
{u,d,r,p,q},
u=N@Exp[σ Sqrt[T/n]];
d=(1/u);
r=Exp[rf T/n];
p=((r-d)/((u-d)r));
q=(1/r)-p;
Sum[Max[S0 u^j d^(n-j)-K,0] Binomial[n,j]p^j q^(n-j),{j,0,n}]
]

binomialPutPrice[S0_,K_,σ_,rf_,T_,n_]:=Module[
{u,d,r,p,q},
u=N@Exp[σ Sqrt[T/n]];
d=(1/u);
r=Exp[rf T/n];
p=((r-d)/((u-d)r));
q=(1/r)-p;
Sum[Max[ K-S0 u^j d^(n-j),0] Binomial[n,j]p^j q^(n-j),{j,0,n}]
]

We use these functions to find the prices of call and put options when the current stock price is $50, the strike price is $45, the annual volatility is 40%, the risk-free rate is 10%, the time to maturity is half a year and the total number of up and down moves is 500.

binomialCallPrice[50,45,0.4,0.1,0.5,500]

binomialPutPrice[50,45,0.4,0.1,0.5,500]

12.2 Black–Scholes–Morton Options Pricing Model

Similar to the binomial option pricing formula defined in the last section, we follow Hull [] to explain the Black–Scholes–Morton option pricing formulas. Define the variables:

(Subscript[c, 0]): the current call option value

(Subscript[p, 0]): the current put option value

(Subscript[S, 0]): the current stock price

K: the exercise price

T: the time in years until the option expires

rf: the annual risk-free interest rate

σ: the annual standard deviation of the rate of return of the underlying stock

δ: the annual dividend yield of the underlying stock

Furthermore, assume that ϕ(x)=((e^(-((x^2)/2)))/(Sqrt[2 π])) is the standard normal density function (where e≈2.71828… the base of the natural logarithms) and let N(x)=(∫_(-∞)^∞)ϕ(x)d x be the standard normal cumulative distribution function, so that N(d) denotes the probability that a random variable drawn from a standard normal distribution is less than d. Then the call and put values can be computed as

(Subscript[c, 0])=(Subscript[S, 0])(e^(-δ T))N((Subscript[d, 1]))-K (e^(-rf T))N((Subscript[d, 2])),

(Subscript[p, 0])=K (e^(-rf T))(1-N((Subscript[d, 2])))-(Subscript[S, 0])(e^(-δ T))(1-N((Subscript[d, 1]))),

where

(Subscript[d, 1])=ln(((Subscript[S, 0])/K))+(rf-δ+((σ^2)/2))T,

(Subscript[d, 2])=(Subscript[d, 1])-σ(Sqrt[T]).

Table  summarizes the price sensitivity measures of call and put options (denoted by Greek symbols) with respect to their major price determinants; here V stands for value of the option.

The built-in Mathematica function FinancialDerivative computes the values and other price sensitivity measures for common types of derivative contracts. The function can compute the value of an option, any of delta, gamma, theta and vega, as well as the implied volatility of the contract. The Mathematica documentation provides many examples of how to use FinancialDerivative. Here are the first 10 of a list of 101 available contracts.

Take[FinancialDerivative[],10]//Column

For example, we compute the price and Greeks of the European-style put option with strike price $50, expiration date 0.3846 years, interest rate 5%, annual volatility 20%, no annual dividend and current price $49.

FinancialDerivative[{"European","Call"},{"StrikePrice"->50,"Expiration"->0.3846},{"InterestRate"->0.05,"Volatility"->0.2,"CurrentPrice"->49,"Dividend"->0},{"Value","Greeks"}]

Similarly, we compute the implied volatility of an American-style call option with the same values of the parameters.

FinancialDerivative[{"American","Call"},{"StrikePrice"->50,"Expiration"->0.3846,"Value"->45},{"InterestRate"->0.05,"CurrentPrice"->49,"Dividend"->0},"ImpliedVolatility"]

One interesting application of FinancialDerivative is to get real-world data and compute related measures related to options. We define computeOptionParameters that computes the theoretical value of options and their Greeks, which takes five arguments:

• ticker symbol (symbol)

• strike price (K)

• expiration date (expirationDate)

• risk-free rate (rf)

• either "European" or "American" (type)

computeOptionParameters[symbol_String,K_,expirationDate_,rf_,type_]:=Module[
{startDate,endDate,returns,σ,price, timeToExpiration,δ,callprice,putprice},
startDate=DatePlus[Today,{-1,"Year"}];
endDate=DatePlus[Today,{-1,"Day"}];
returns=QuantityMagnitude[FinancialData[symbol,"FractionalChange",{startDate,endDate}]["Values"]];
σ=StandardDeviation[returns]Sqrt[252]/100;
price=FinancialData[symbol,"Price"][[1]];
timeToExpiration=DateDifference[Today,DateObject[expirationDate]][[1]]/360//N;
δ=Total@QuantityMagnitude[FinancialData[symbol,"Dividend",{startDate,endDate}]["Values"]]/price;
{callprice,putprice}=FinancialDerivative[{type,#},{"StrikePrice"->K,"Expiration"->timeToExpiration},{"InterestRate"->rf,"Volatility"->σ,"CurrentPrice"->price,"Dividend"->If[NumberQ@δ,δ,0]},{"Value","Greeks"}]&/@{"Call","Put"};
ArrayFlatten[
{
{{#}&/@Flatten[{" ","Price",Keys[putprice[[2]]]}],
Prepend[
Transpose[{
Flatten[{callprice[[1]],Values@callprice[[2]] }],
Flatten[{putprice[[1]],Values@putprice[[2]] }]
}],
{"Call","Put"}]}
}
]
]

We compute the European-style option parameters for Boeing (BA), assuming that the option expires on December 28, 2020, with exercise price $145 and risk-free rate 0.0187. Make sure that the expiration date is later than the current date, since computeOptionParameters uses historical data.

Grid[
computeOptionParameters["BA",320,{2020,12,28},0.0187,"European"],
Frame->All]//Text

We strongly encourage you to explore the built-in or online documentation for the powerful FinancialDerivative function.

13. Conclusion

The article provides a brief overview of built-in functions and introduces many functions especially designed for analysis of financial data. In particular, we have focused on functions that are more relevant to introductory computational financial concepts. We emphasize importing company fundamental data and its visualization, analysis of individual stocks and portfolio returns, factor models and the use of built-in functions for bond and financial derivative analysis. The functions we have provided are just a few examples. The Wolfram Language can do much more than what we have shown in this article. Interested readers can start exploring the Wolfram Language via Mathematica’s extensive documentation.

References

H. R. Varian, ed., Economic and Financial Modeling with Mathematica, New York: Springer-Verlag, 1993.

H. R. Varian, ed., Computational Economics and Finance Modeling and Analysis with Mathematica, New York: Springer-Verlag, 1996.

W. Shaw, Modelling Financial Derivatives with Mathematica, Cambridge, UK: Cambridge University Press, 1998.

S. Stojanovic, Computational Financial Mathematics using MATHEMATICA Optimal Trading in Stocks and Options, Boston: Birkhäuser, 2003.

F. Black, “Capital Market Equilibrium with Restricted Borrowing,” Journal of Business, 4(3), 1972 pp. 444–455. www.jstor.org/stable/2351499.

G. Campolieti and R. Makarov, Financial Mathematics: A Comprehensive Treatment, London: Chapman and Hall/CRC Press, 2014.

E. F. Fama and K. R. French, “A Five-Factor Asset Pricing Model,” Journal of Financial Economics, 116(1), 2015 pp.1–22. doi:10.1016/j.jfineco.2014.10.010.

C. R. Bacon, Practical Risk-Adjusted Performance Measurement, 2nd ed., Hoboken: Wiley, 2013.

Fidelity Learning Center. “Technical Indicator Guide.” (Jul 29, 2020) www.fidelity.com/learning-center/trading-investing/technical-analysis/technical-indicator-guide/overview.

S. Benninga, Financial Modeling, 4th ed., Massachusetts: The MIT Press, 2014.

J.C. Hull, Options, Futures and Other Derivatives, 10th ed., New York: Pearson Education Limited, 2018.