Mathematica® Beyond Mathematics: The Wolfram Language in the Real World

11 Economic and Financial Applications

The latest Mathematica versions include multiple capabilities related to economics and finance. The new functionality makes it easy to perform financial visualizations and to access both historical and real time data. The commands related to stock markets and financial derivatives would be especially useful to readers interested in investment portfolio management. Additionally, the functions for optimization, with applications in economics and many other fields, have been improved significantly. A new function to solve the Traveling Salesman Problem in a very efficient way has been included. All of these topics will be discussed in the sections below using real world examples. This chapter should be read along with Chapter 6, the one covering probability and statistics.

11.1 Financial Information

11.1.1 Introduction

Mathematica, as part of its computable data functionality, includes the command FinancialData to access finance-related information. Although in theory the program can display thousands of indicators, in reality the accessible ones are mostly those related to the US markets. Those readers who may need reliable access to other markets or real time financial data would probably be interested in taking a look at the Wolfram Finance Platform, a new product from Wolfram Research (the developers of Mathematica) that focuses on finance: http://www.wolfram.com/finance-platform/.

In any case, although the information may not be directly accessible from FinancialData, nowadays it would be easy to access it. For example, many bank platforms enable their clients to download data in an Excel format that could then be imported into Mathematica and analyzed using any of the multiple functions related to finance:

http://reference.wolfram.com/mathematica/guide/Finance.html

11.1.2 FinancialData

FinancialData ["name", "property", {start, end, interval}] gives the value of the specified property for the financial entity “name” (it could be an index, a stock, a commodity, etc.) for a given period of time (by default it returns the most recent value). We recommend to read the documentation pages to know more about this function’s capabilities.

The symbols (“name”) used are the ones available in http://finance.yahoo.com/ since Yahoo Finance supplies the information to the function. In the website we can find that the symbol for the Standard & Poor’s 500 index, that includes the 500 largest public companies in terms of market capitalization that trade in either NYSE or NASDAQ is: ^GSPC (we can also use SP500). We can check it as follows:

FinancialData["SP500", "Name" ]

S&P 500
If we want to see the price (with a few minutes delay) we can type:

FinancialData["^GSPC"]

2269.
We can also access historical data. In this example we show the trajectory of the S&P 500 index from January 1, 2000 to December 31, 2016.

DateListPlot[FinancialData ["^GSPC",
{"January 1 2000", "December 31 2016" }], Joined → True]

If the financial entity is trading at the time of the query, the price that we get is the one from a few minutes earlier (usually the delay is 15 minutes). The access to information in “real time” requires additional payments to 3rd party providers such as Bloomberg or Reuters. There are even people willing to pay very large amounts of money to gain an advantage of a few milliseconds when trading:

http://adtmag.com/articles/2011/07/29/why-hft-programmers-earn-top-salaries.aspx.

The following function will dynamically update the S&P500 index (as long as the command is executed during trading hours):

sp500 : = FinancialData["SP500", "LatestTrade"]
The output changes every two seconds:

Dynamic[sp500, UpdateInterval → 2]

sp500
We can also represent it in a self-updating graph:

data = {};
Dynamic[Last[AppendTo[data, sp500]],
UpdateInterval → 2, TrackedSymbols → {}]

Dynamic[DateListPlot[data, Joined → True], SynchronousUpdating → False]
As mentioned earlier, if the entity is not available in FinancialData we may still be able to import the information from specialized websites. A very interesting one is Quandl, http://www.quandl.com, that provides access to numerous economic and financial databases. Most of them are available for free (https://www.quandl.com/search?type=free) but some of them require subscription (https://www.quandl.com/search?type=premium). We can first select the data and download it in our chosen format: Excel, CSV or any other. Then we can use Import to analyze them using Mathematica. Alternatively, we can access the data directly using Import or URLFetch. If you use Quandl frequently, you may be interested in downloading QuandlLink, a free Mathematica package containing the function QuandlFinancialData, similar to FinancialData. The package is available in: https://github.com/bajracha71/Quandl-Mathematica-QuandlLink (Accessed on December 4, 2016).

Another interesting feature is the ability to generate email messages linked to certain events. For further information please visit: guide/AutomaticMailProcessing.

For example, we may be interested in sending a message if a stock price goes above or below a certain threshold.

First we need to provide Mathematica with our email account information. This is done in Edit ▶ Preferences... ▶ Internet Connectivity ▶ Mail Settings. In the case of Hotmail, the configuration is the one shown in Figure 11.1. For Gmail is similar with “Server” → “smtp.gmail.com” and “PortNumber” → 587. For other popular email services such as Yahoo, you can find the relevant information on the Internet. For enterprise mail servers you may have to ask your system administrator.

Figure 11.1 Mail Settings in Mathematica.
Now we can send a test message:

SendMail["To" -> "myemail@gmail.com",
"Subject" -> "Test message",
"Body" -> "Test message",
"Password" -> "your password"]
If it worked correctly, we can now define the following function:

ss[mens_] := SendMail["To" -> "myemail@hotmail.com",
"Subject" -> "Message via sendmail",
"Body" -> ToString[mens],
"Password" -> "your password"]
You can use this function at your convenience. In this example, we create a stopwatch and send a message when it reaches 20. You can create similar commands to send messages when a stock goes above or below a certain threshold.

Dynamic[Refresh[If[Round[Clock[1000]] == 20,
ss[Clock[1000]], Round[Clock[1000]]], UpdateInterval → 1]]
In the next example we send a message at a scheduled time. The command below activates the task, the latest Euro to US Dollar exchange rate, 10 seconds from the current time, and sends an email 1 second after that.

RunScheduledTask[ss[FinancialData["EUR/USD"]],
{1}, AbsoluteTime[DateString[]] + 10]

11.1.3 Visualization

The number of specialized graphs for economic and financial visualization has increased significantly in the latest versions of Mathematica.

All the functions below end in Chart. We have already seen some of them in Chapter 6. In this section, we’re going to cover the ones normally used for analyzing stock prices.

?*Chart

▼ System

Choose any of them and click on it to access its documentation.

These graphs complement the capabilities of FinancialData but they can also be used with other functions. Let’s see some of them:

Let’s visualize the historical closing prices for the S&P 500 using different graphs.

data = FinancialData["SP500", "Close", {{2016, 1, 1}, {2016, 07, 31}}];
If we click on any point inside the plot, we’ll see the information associated to it.

GraphicsRow[
{KagiChart[data], PointFigureChart[data], LineBreakChart[data]}]
The next graph uses the function RenkoChart with EventLabels and Placed to add symbols (“→” y “↓”) indicating in what position they must be placed.

RenkoChart[{"NYSE:GE", {{2016, 01, 1}, {2016, 07, 31}}}, EventLabels →
{{2015, 10, 1} → Placed["→", Before], {2016, 5, 2} → Placed["↓", Above]}]
The graph below displays the trajectory of GE in the NYSE from January 1 to July 31, 2016. For each trading day, we can see the open, high, low and close prices.

CandlestickChart[{"NYSE:GE", {{2016, 01, 1}, {2016, 7, 31}}}]
ChartElementFunction is an option that when added to BarChart or PieChart3D, improves the presentation of the data.

GraphicsRow[{BarChart[{{1, 5, 2}, {4, 5, 7}},
   ChartElementFunction → "GradientScaleRectangle"],
  PieChart3D[{3, 2, 2, 4, 1}, ChartElementFunction → "TorusSector3D",
   SectorOrigin → {Automatic, 1}, ChartStyle → "Rainbow"]}, Frame → All]

In the world of professional stock market investors, there’s a group that uses chart analysis as its main tool for making investment decisions. Its members are known as chartists and their techniques are known as technical analysis. They closely study the past performance of stocks or indices to find trends and support and resistance levels. These people associate certain graph shapes to future markets ups and downs. Although there’s no sound mathematical theory justifying their predictions, all economic newspapers have a section dedicated to this type of analysis.

The powerful function InteractiveTradingChart, introduced in Mathematica 8, includes practically all the functionality that a chartist may require. In this example, we display GE’s share price history for a given period. Notice that in the upper area of the graph, we can see a trimester at a time and using the bottom part we can move through the entire period. We can also choose the time interval (days, weeks or months) and even what chart type and indicators to show.

InteractiveTradingChart[{"NYSE:GE", {{2016, 01, 1}, {2016, 07, 31}}}]

11.1.4 Automatic Adjustments

With Mathematica we can automatically fit our financial data to a given distribution (by default the program will use the maximum likelihood method) with the function EstimatedDistribution.

Let’s fit Google’s share price to a lognormal distribution after downloading the data.

googleStock =
FinancialData["GOOG", {{2010, 3, 10}, {2016, 07, 31}, "Day"}, "Value"];

dist = EstimatedDistribution[googleStock, LogNormalDistribution[μ, σ]]

LogNormalDistribution[6.03238, 0.360796]
The fit is poor on both tails. As a matter of fact, all the efforts that have been made trying to predict financial markets have produced unsatisfactory results. These markets exhibit what Mandelbrot (the father of fractal geometry) called “wild randomness” (single observations can impact the total in a very disproportionate way).

QuantilePlot[googleStock, dist, Filling → Automatic, FillingStyle → Red]

11.2 Financial Functions

The functions related to financial calculations can be found in guide/Finance. We will give an overview of some of them in this section.

11.2.1 Time Value Computations

Clear["Global`*"]

The following functions are included under the Time value Computations category: TimeValue, EffectiveInterest, Annuity, AnnuityDue and Cashflow. Let’s see some examples.

TimeValue ["s", "i", "t"] calculates the time value of a security s at time t for an interest specified by i.

If we invest an amount c during t years at an annual interest rate of i, to find the future value Mathematica will apply the following formula:

TimeValue[c, i, t]

c (1 + i)^t
How much do we need to invest at 3.5% annual interest rate to get €100,000 in 10 years?

Solve[TimeValue[quantity, .035, 10] == 100 000, quantity]

{{quantity → 70 891.9}}
We have two options to invest €60,000 for 10 years. Option A offers a fixed annual interest rate of 4% during the 10 years. Option B offers an interest rate of 3% during years 1 and 2, 4 % in years 3 to 5, and 4.5% until the end of the investment period. Which option is the best one?

{TimeValue[60 000, 0.04, 10] (*product A*),
TimeValue[60 000, {{0, 0.03}, {2, 0.04}, {5, 0.045}}, 10](*product B*)}

{88 814.7, 89 229.2}

The function EffectiveInterest ["r", "q"] returns the effective interest rate corresponding to an interest specification r, compounded at time intervals q. It is very useful to avoid making mistakes when comparing nominal and effective interest rates. Annuity ["p", "t"] represents an annuity of fixed payments p made over t periods. Below we show examples using both functions.

What is the effective rate of an annual nominal interest rate of 4% compounded quarterly?

EffectiveInterest[.04, 1/4]

0.040604
We’d like to deposit €60,000 in an account with an annual interest of 4% compounded quarterly. How much money will we be able to withdraw every quarter so the money lasts for 10 years?

Solve[TimeValue[Annuity[{wd, {-60 000}}, 10, 1/4],
EffectiveInterest[.04, 1/4], 10] == 0, wd]

{{wd → 1827.34}}
The next set of functions creates a loan amortization table for a given period and interest rate. In this case we borrow €150,062 and make monthly payments of €1,666 during 10 years at 6% nominal annual rate. The same approach can be used with different numbers:

years = 10; monthlypayments = 1666; monthlyinterest = 0.06/12;

loanValue[month_] = TimeValue[
   Annuity[monthlypayments, years*12 - month], monthlyinterest, 0];

amortizationList = Table[
    {i + 1, loanValue[i], monthlypayments, loanValue[i] *monthlyinterest,
     monthlypayments - loanValue[i] *monthlyinterest,
     loanValue[i + 1]}, {i, 0, years*12 - 1}] // Round;

amortizationTable[month1_, month2_] :=
TableForm[amortizationList[[month1 ;; month2]], TableHeadings →
   {None, {"Month", "Beginning\nPrincipal", "Monthly\nPayment",
     "Interest\nPayment", "Principal\nPayment", "Ending\nPrincipal"}}]
Let’s see the amortization table during the first 12 months. Try to display the entire period of the loan.

amortizationTable[1, 12]
We can check graphically how the remaining principal goes down over time until it becomes 0 at the end of the loan:

ListLinePlot[Table[{i + 1, loanValue[i]}, {i, 0, years*12 - 1}],
AxesLabel → {"Month", "Remaining principal"}, ImageSize → 300]

11.2.2 Bonds

Clear["Global`*"]

Bonds are financial assets that we can also analyze with Mathematica. When an organization needs capital, instead of borrowing from banks, they can issue bonds. Bonds can be considered a type of loan in which the lenders are the ones purchasing them. Those lenders can usually be anyone. There are different types of bonds. Normally they pay a fixed amount of interest at predetermined time intervals: annually, half-annually, quarterly or even monthly. Some bonds issued for short periods of time (less than one year), such as Treasury Bills, don’t pay interest explicitly, they are sold at a discount from their par value.

For analyzing bonds, Mathematica has the function FinancialBond. It takes the following arguments (some of them are optional):

“FaceValue”	face value, par value
“Coupon”	coupon rate, payment function
“Maturity”	maturity or call/put date
“CouponInterval”	coupon payment interval
“RedemptionValue”	redemption value
“InterestRate”	yield to maturity or yield rate
“Settlement”	settlement date
“DayCountBasis”	day count convention

Let’s see some examples.

The yield to maturity of a €1,000 30-year bond maturing on June 31, 2018, with a coupon of 6% paid quarterly that was purchased on January 6, 2012 for €900 would be:

FindRoot[FinancialBond]{"FaceValue" → 1000, "Coupon" → 0.05,

"Maturity" → {2018, 6, 31}, "CouponInterval" → 14
{"InterestRate" → y, "Settlement" → {2012, 1, 6}}] == 900, {y, .1}]

{y → 0.069269}
A zero-coupon bond with a redemption value of €1,000 after 10 years is currently sold for €400. Find the implied yield rate compounded semiannually:

FindRoot[FinancialBond[

{"FaceValue" → 1000, "Coupon" → 0, "Maturity" → 10, "CouponInterval" → 12},

{"InterestRate" → y, "Settlement" → 0}] == 400, {y, .1}]

{y → 0.0937605}

11.2.3 Derivatives

Clear["Global`*"]

In Chapter 41 of the book of Genesis, the Egyptian Pharaoh asks Joseph to interpret a dream he’s had about 7 sleek and fat cows followed by 7 ugly and gaunt ones. Joseph interprets the dream as periods of plenty followed by periods of famine and gives the Pharaoh the following piece of advice: “... Let Pharaoh appoint commissioners over the land to take a fifth of the harvest of Egypt during the seven years of abundance. They should collect all the food of these good years that are coming and store up the grain under the authority of Pharaoh, to be kept in the cities for food. This food should be held in reserve for the country, to be used during the seven years of famine that will come upon Egypt, so that the country may not be ruined by the famine”. Some people consider this story as the earliest mention of the derivative concept based on the idea of fixing in advance the future price of a certain item. Closer to our current era, on April 3, 1848, the Chicago Board of Trade was established to create a marketplace for the buying and selling of several commodities (grains, cattle, hogs, etc.) at fixed prices and following certain standards. Normally, the buyer would pay in advance a small portion of the total amount to the seller. The purpose was clear: the buyer could set in advance how much things were going to cost in the future and the seller could secure a certain selling price; it could be considered as an insurance against uncertainty. Since then, this idea has been applied to many other areas: raw materials, shares, stock indexes, etc.

In finance, a derivative product is one whose price depends on the price of another product (the underlying asset). We say that the price of the former derives from the price of the latter. For example, we can create a derivative based on the price of gold without the need to physically own the gold. In principle, these instruments should play a beneficial economic role. Very often, however, people associate them with speculative transactions more similar to gambling than to uncertainty reduction. There’s a kernel of truth is this opinion since investments in these products can generate either large losses or large gains. This is the reason why there have always been investors trying to (and sometimes succeeding) influence the price of the underlying for speculative purposes. Warren Buffett, the most successful investor of the 20th century and one of world’s richest men, even described them as “weapons of mass destruction” although it’s known that some of his investments take advantage of them, and that if used properly they help to reduce uncertainty.

When considering derivatives as an alternative to purchasing stocks, the two most important things to keep in mind are: i) For the same number of shares, derivatives require a much smaller investment. As a consequence, gains and losses are significantly higher, ii) When the share price goes up, everybody wins (more precisely, the latent value is the same for everybody). With derivatives, when someone wins somebody else loses and the other way around (it’s a zero-sum game).

An example of a financial derivative is a stock option, or a well-known variant of it more commonly available to regular investors, a warrant. An option or warrant gives the owner the right to buy (Call) or sell (Put) the underlying asset at a given price usually determined at the beginning of the contract. Both have prices (premiums) that are normally much smaller than the stock associated to them. For example, we could buy a warrant to purchase Telefónica shares in a year for €15 each and pay €0.5 for it. At expiration, if the shares were trading above €15, let’s say €16.5, we would make a profit of €1 per share, after subtracting the €0.5 initially paid for it (it’s common to just settle the differences, without actual delivery of the underlying asset). However, if at expiration Telefónica was trading below €15, we would just lose the premium paid.

An informative overview of financial derivatives and their statistical measurement can be found in:

https://www.imf.org/external/pubs/ft/wp/wp9824.pdf.

The list below contains some of the most commonly used terms when describing derivatives transactions:

Underlying asset: the financial instrument on which a derivative’s price is based.

Call option. An agreement that gives the buyer the right (but not the obligation) to buy the underlying asset at a certain time and for a certain price from the seller.

Put option. An agreement that gives the buyer the right (but not the obligation) to sell the underlying asset at a certain time and for a certain price to the seller.

Premium. The amount paid for the right to buy or sell the underlying asset.

Market price. The current price at which the underlying asset can be bought or sold.

Strike. The price at which the agreement can be exercised.

The Greek letters delta, gamma, theta and vega are used to describe the risk of changes in the underlying asset.

Delta: It’s the option sensitivity to changes in the underlying asset. It indicates how the price of the option changes when the underlying price goes up or down. For an option at the money (strike price same as market price), delta is usually around 50%. This means that if the price of the underlying asset changes by €2, the price of the option will change by €1.

Gamma: It measures the sensitivity of delta to price changes in the underlying asset. Since delta is not a linear function, gamma is not linear either. Gamma decreases as the uncertainty increases.

Vega: This metric indicates how sensitive the option price is to volatility changes in the underlying asset. It represents the change in the option price if the volatility of the underlying asset changes by 1%. When volatility goes up, the price of the option also goes up. Vega may vary depending on whether the option is at the money, in the money or out of the money.

Theta: A measure of the decrease in the value of an option as a result of the passage of time (easy to remember since both start with t). The longer the time period until exercise, the higher Theta is: There’s more time for the underlying asset’s price to move in a direction favorable for the option holder.

Mathematica gives us a very wide array of options when using the function:

FinancialDerivative

["instrument", "parameters", "ambient parameters", "property (optional)"] returns the value of the specified financial instrument, computed for the stated property.

To make the valuation examples in this section as realistic as possible we will need certain market data.

Often, to make the purchase of an option cheaper, we can sell another one in such a way that we limit our potential gains but we also reduce the premium paid. This is the idea behind zero cost strategies.

An example of this strategy would be as follows: We purchase an option with a strike price equal to the exercise price (an at the money option) and sell another one with a higher strike price, let’s suppose 18% higher. Then, this figure will be our maximum return. Obviously, the higher the potential return, the higher the potential risk and the other way around. This strategy is known as a Call Spread.

In this case we use Telefónica shares trading in the NYSE (the data and functions used in the example come from Sebastián Miranda from Solventis). The stock prices for the Spanish stock exchange (IBEX) are currently not available.

FinancialData["NYSE:TEF", "Name"]

Telefonica SA
We’re going to need the 12-month EURIBOR rate. This type of interest rate is published by Euribor EBF (http://euribor-rates.eu). Let’s import the 10 most recent rates (remember that in Chapter 2 we explained how to extract specific parts of a web page) as follows:

Import["http://www.euribor-rates.eu/euribor-rate-12-months.asp", "Data"][[
2, 4, 1, 2, 3]]

{Rate on first day of the year,
{{01-04-2016, 0.058%}, {01-02-2015, 0.323%},
  {01-02-2014, 0.555%}, {01-02-2013, 0.543%}, {01-02-2012, 1.937%},
  {01-03-2011, 1.504%}, {01-04-2010, 1.251%}, {01-02-2009, 3.025%},
  {01-02-2008, 4.733%}, {01-02-2007, 4.030%}}}
However, unless we’re creating a program to automate the process, probably the easiest thing to do is to copy the values directly. Besides, it’s unusual for websites to keep their structure for a long time so, by the time you read these lines, the instruction above may not work. We’ll store in memory the first rate for 2016 that we will use for our calculations later on.

euribor12m = 0.058/100

0.00058
For our computations, we need to know several indicators associated to the stock. The values returned are time sensitive and depend on the time of execution of the query. We show the commands below so readers can replicate the example. Notice that after execution, the predictive interface gives us the possibility of making other queries (Figure 11.2).

FinancialData["NYSE:TEF", "Price"] (* stock price,
in USD, in NYSE at the time of the query*)

9.87

Figure 11.2 Predictive interface for financial data.

spotPrice = 9.87;
Stock volatility during the past 50 days. (Volatility is a measure of the variability of the stock price over a specific time horizon. It’s usually the standard deviation of the closing stock price during the last 50 trading sessions, or an alternative chosen time horizon):

FinancialData["NYSE:TEF", "Volatility50Day"]

0.528374

vol = 0.528374;
Last year’s dividend yield, computed as the ratio between the dividends paid and the current stock price.

FinancialData["NYSE:TEF", "DividendYield"]

0.0904

div = 0.0904;
Now we can evaluate the price of a call option given its expiration date (in this case, in one year), and using the current stock price as the strike price:

FinancialDerivative[{"American", "Call"},
{"StrikePrice" → spotPrice, "Expiration" → 1},
{"InterestRate" → euribor12m, "Volatility" → vol,
"CurrentPrice" → spotPrice, "Dividend" → div}]

1.67404
In the next graph, we can see the evolution of the value of our call spread as a function of the underlying asset price and the remaining time. To make it easier to visualize, the solid green line indicates the current price, our strike price; and the dashed red line, the price of the underlying asset at which we maximize our profit. We need to keep in mind that our starting point in the time scale is 1. Remember that our goal is to have a maximum gain of 18%.

data =
  {#, 1- #2, FinancialDerivative[{"American", "Call"}, {"StrikePrice" →
         spotPrice, "Expiration" → 1}, {"InterestRate" → euribor12m,
        "Volatility" → 0.5, "CurrentPrice" → #, "Dividend" → div,
        "ReferenceTime" → #2}] - FinancialDerivative[{"American", "Call"},
       {"StrikePrice" → 1.18*spotPrice, "Expiration" → 1},
       {"InterestRate" → euribor12m, "Volatility" → 0.5,
        "CurrentPrice" → #, "Dividend" → div, "ReferenceTime" → #2}]} & @@@
   Flatten[Table[{price, referencetime}, {referencetime, 0, 1, 0.05},
     {price, 0.8 spotPrice, 1.4*spotPrice, 1}], 1];

ListPlot3D[data,
(* Contours *)
Mesh → {{{spotPrice, {Green, Thick}},
    {1.18*spotPrice, {Red, Thick, Dashed}}}},
(* Formats *)
AxesLabel →
   (Style[#, 12] & /@ {"Underlying Price", "Remaining Time", "Profit"}),
ColorFunction → ColorData["LightTemperatureMap"],
ColorFunctionScaling → True]
At expiration, the profile of our profit as a function of the underlying asset price is:

ListPlot[Cases[data, {x_, 0., y_} → {x, y}],
Joined → True, PlotStyle → Directive[Thickness[0.02], Red],
AxesLabel → {"Underlying Price", "Profit"}]

In this example, the warrants that we’ve used are of the American type, meaning that they can be exercised at any time until the expiration date, although it is not usually optimum to do so.

However, Mathematica, in addition to financial derivatives price calculations, can also compute their sensitivities, known as Greeks. These sensitivities are the partial derivatives of the option price with respect to different valuation parameters. For example, the delta of an option is the partial derivative of its price with respect to the price of the underlying asset. This measure not only gives us an idea of how sensitive the option value is to small changes in the price of the underlying asset but also, for standard options (also known as vanilla), it gives us the probability that the option will be exercised. Gamma, the second derivative with respect to the underlying asset price, tells us how sensitive the delta is to changes to it. This metric is very important for hedging purposes. Rho is the partial derivative with respect to the interest rate and in the case of stock options is the least sensitive Greek. Theta is the derivative with respect to time to expiration. Initially, since options are not usually bought in the money i.e. if Telefónica trades at €12, purchasing an option for €10 is not very common, they don’t have intrinsic value. They only get value from their time component. The final one, Vega, is the partial derivative with respect to the volatility of the underlying asset. This is such an important parameter in option valuation that for certain options, their prices are quoted in terms of the volatility of their respective underlying assets.

The graph below shows how the Greeks change with respect to the remaining time to expiration:

ticks[min_, max_] := Table[
  {j, ToString[Round[12 - j 12]] <> " Months", {.03, 0}}, {j, min, max, 0.2}]

TabView[Rule @@@ Transpose[{{"Delta", "Gamma", "Rho", "Theta", "Vega"},
    ListPlot[#, Ticks → {ticks, Automatic}] & /@
     Transpose[Table[{i, #} & /@ FinancialDerivative[{"American", "Call"},
          {"StrikePrice" → 55.00, "Expiration" → 1},
          {"InterestRate" → 0.1, "Volatility" → 0.5, "CurrentPrice" → 50,
           "Dividend" → 0.05, "ReferenceTime" → i},
          "Greeks"][[All, 2]], {i, 0, 0.95, 0.05}]]}]]

The same Greek letters are used to indicate the risks of portfolio of options due to changes in the price of the underlying assets or market conditions.

Delta: the delta of a portfolio of options with the same underlying assets is the weighted average of the delta of the individual options.

The delta factor = the change in the price of the option as a result of changes in the price of the underlying assets. For an at-the-money option, the delta factor is usually close to 50%, so if the underlying price changes by €1, the price of the option will change by €0.5.

11.2.4 The Black–Scholes Equation

In this example we first compute the price of a derivative using Mathematica’s existing function and then compare it with the result obtained by solving the Black–Scholes equation directly

(http://www.wolfram.com/language/11/partial-differential-equations/find-the-value-of-a-european-call-option.html).

Let’s start by using FinancialDerivative to compute the price of a European vanilla call option with the data from the previous section:

spotPrice = 9.87; euribor12m = 0.00058; vol = 0.528374;
FinancialDerivative[{"European", "Call"}, {"StrikePrice" -> spotPrice,
"Expiration" -> 1}, {"InterestRate" -> euribor12m, "Volatility" -> vol,
"CurrentPrice" -> spotPrice}]

2.05882
Now, let’s do the same computation using the Black–Scholes equation (To create more visually appealing documents you can convert the input cell to the traditional format: ▶ Convert to ▶ TraditionalForm):

BlackScholesModel= {−r c(t,s)+r s∂ c(t,s)∂s+12s2σ2∂2 c(t,s)∂s2+∂ c(t,s)∂t=0, c(T, s)=max(s−k,0)};
After solving the boundary value problem we get:

(dsol = c[t, s] /.
DSolve[BlackScholesModel, c[t, s], {t, s}][[
1]]) // TraditionalForm

12e−rT(S erT erfc((2 log(k)+(2r+σ2)(t−T)−2 log(s))/(22 σT−t))−k ert erfc((2 log(k)+(2r−σ2)(t−T)−2 log(s))/(22 σT−t)))
Finally, we calculate the price of the derivative and compare it with the one previously obtained. We can see it’s the same:

dsol /. {t -> 0, s → spotPrice, k -> spotPrice, σ -> vol, T -> 1,
r -> euribor12m}

2.05882

11.2.5. Basel II Capital Adequacy: Internal Ratings-Based (IRB) Approach

The key cornerstone of prudential regulation of banks (and financial institutions in general) is to ensure that each bank holds sufficient equity capital to absorb unexpected losses, that is, the materialization of financial risks, mainly market (price), credit, and operational risks.

The example below, Basel II Capital Adequacy: Internal Ratings-Based (IRB) Approach, by Poomjai Nacaskul shows how to calculate, using the Internal Ratings-Based approach, the minimum amount of capital that a bank should keep following the standards defined by the Basel II regulatory framework on bank capital adequacy.

http://demonstrations.wolfram.com/BaselIICapitalAdequacyInternalRatingsBasedIRBApproach

To facilitate the understanding of the example, we also include the code:

(*Setting defaults*)
DefaultEAD = 1 000 000;
MinimumPD = 0.0001; MaximumPD = 0.25;
MinimumLGD = 0.01; DefaultLGD = 0.7; MaximumLGD = 1.0;
MinimumMaturity = 1/250;
DefaultMaturity = 1;
MaximumMaturity = 30;

(* Create shorthands for Normal CDF & CDF inverse *)
Phi[x_] := CDF[NormalDistribution[0, 1], x];
PhiInverse[p_] := InverseCDF[NormalDistribution[0, 1], p];
(* Correlation Function, Basel Committee's way of parameterising a 2-
parameter formula (PD, rho) with 1 parameter (PD). *)
R[PD_] := 0.12* (1 - Exp[-50 PD]) / (1 - Exp[-50]) +
   0.24* (1 - (1 - Exp[-50 PD]) / (1 - Exp[-50]));
ConditionalPercentageExpectedLoss[PD_, LGD_] := Phi[
   (PhiInverse[PD] + Sqrt[R[PD]] PhiInverse[0.999]) /Sqrt[1 - R[PD]]] *LGD
UnconditionalPercentageExpectedLoss[PD_, LGD_] := PD*LGD;

(* Maturity Adjustment, with a 1-year maturity as a base case *)
b[PD_] := (0.11852 - 0.05478*Log[PD])^2;
MaturityAdjustment[Maturity_, PD_] :=
  (1 + (Maturity - 2.5) *b[PD]) / (1 - 1.5*b[PD]);

(* Capital Requirement Calculation: note as PD --> 1,
this quantity --> 0, as credit risk is taken care of via provisioning,
i.e. E[Loss] = exposure x PD x LGD. *)
CapitalRequirement[PD_, LGD_, Maturity_] :=
  (ConditionalPercentageExpectedLoss[PD, LGD] -
     UnconditionalPercentageExpectedLoss[PD, LGD]) *
  MaturityAdjustment[Maturity, PD];

Manipulate[
Pane[Column[{Text@Style["Basel II Capital Adequacy", Bold, 16],
    Text@"Regulatory (minimal) capital for financial institutions
       per $ exposure,\nas per Internal Ratings-Based
       Approach (IRB) Risk Weight Function (RWF)\n",
    Grid[{{Plot[CapitalRequirement[PD, LGD, Maturity], {PD, 0.0001, 0.25},
        PlotLabel -> "as a function of Probability of Default (PD)\nat
           a given Maturity (M)", ImageSize → 320],
       Plot3D[CapitalRequirement[PD, LGD, Maturity], {PD, 0.0001, 0.25},
        {Maturity, 1/250, 30}, PlotLabel → "as a function of PD and M\n",
        ImageSize → 280, Mesh → None, ColorFunction →
         Function[{x, y, z}, ColorData["GrayTones"][z]]]}}]},
    Alignment → Center], ImageSize → {610, 350}],
{{LGD, 0.7, "Loss Given Default (LGD)"}, 0.01, 1, Appearance → "Labeled"},
{{Maturity, 1, "Maturity M in years (left graph only)"},
  {1/250, 1/50, 1/12, 1/4, 1, 2, 5, 10, 15, 30},
  ControlType → RadioButton}, SaveDefinitions → True]

images

To learn more about financial modeling and risk management using Mathematica, Igor Hlivka has created several excellent example notebooks. You can find them in:

http://community.wolfram.com/web/ihcomm.

11.3 Optimization

11.3.1 What is Constrained Optimization?

The optimization problems that we discuss in this section are related to finding an optimum (maximum or minimum, depending on the case) for an objective function of n variables, given r constraints.

This type of problems is common in several fields such as economics, finance and engineering as we will show using different examples. You can find an excellent tutorial in the documentation pages: “Introduction to Constrained Optimization in the Wolfram Language” (tutorial/ConstrainedOptimizationIntroduction).

To solve them, Mathematica has the following functions:

Function	Application	Description
FindMinimum/FindMaximum	Local Numerical optimization	Linear programming methods, non- linear, interior points, and the use of second derivatives
NMinimize / NMaximize	Global Numerical optimization	Linear programming methods, Nelder- Mead, differential evolution, simulated annealing, random search
Minimize / Maximize	Exact Global optimization	Linear programming methods, algebraic cylindrical decomposition, Lagrange multipliers and other analytical methods, integer linear programming.
LinearProgramming	Linear optimization	Simplex, modified simplex, interior point

The following tree can help us decide the most appropriate function to find the minimum (To find the maximum we would use the functions Maximize, NMaximize or FindMaximum instead).

The command below draws the tree. When creating a document, you may be interested in seeing only the result. To hide the input just click directly on the output marker cell.

TreePlot[{{"Is your problem linear?" → "LinearProgramming", "Yes"},
   {"Is your problem linear?" → "Do you want a global optimum?", "No"},
   {"Do you want a global optimum?" → "Do you want an exact solution?",
    "Yes"}, {"Do you want a global optimum?" → "Is your problem small?",
    "No"}, {"Is your problem small?" → "NMinimize", "Yes"},
   {"Is your problem small?" → "FindMinimum", "No"},
   {"Do you want an exact solution?" → "Minimize", "Yes"},
   {"Do you want an exact solution?" → "NMinimize", "No"}},
Automatic, "Is your problem linear?", DirectedEdges -> True,
  VertexLabeling → True, VertexCoordinateRules ->
   {{0.6358, 2.049}, {-0.26, 1.3}, {1.06, 1.3}, {0.448, 0.54},
    {1.78, 0.54}, {0.76, -0.208}, {1.66, -0.208}, {0, -0.208}}] /.
((# → Button[TextCell[#], Inherited, BaseStyle → Dynamic[If[
            CurrentValue["MouseOver"], {"Link", FontColor →
              RGBColor[0.854902, 0.396078, 0.145098]}, {"Link"}]],
        ButtonData -> "paclet:ref/" <> #, Appearance → "None"]) & /@
  {"LinearProgramming", "NMinimize", "FindMinimum", "Minimize"})

images

In many cases the problem can be solved using several functions. With the exception of LinearProgramming, the syntax for the rest of the functions is very similar: f[{objective function, constraints}, {variables}].

Clear["`Global`*"]

Let’s see some examples.

Given the objective function (ob) with constraints c1 and c2:

var = {x, y}, ob = x + y, c1 = 0 ≤ x ≤ 1, c2 = 0 ≤ y ≤ 2
The maximum can be found using any of the previously mentioned commands:

{NMaximize[{ob,c1,c2}, var], Maximize[{ob,c1,c2}, var],
FindMaximum[{ob,c1,c2}, var]}

{{3., {x → 1., y → 2.}}, {3, {x → 1, y → 2}}, {3., {x → 1., y → 2.}}}
We can see that the maximum occurs when {x → 1, y → 2} and is 3. We can verify it (remember that “/.” is used for substitutions).

x + y /. {x → 1, y → 2}

3
The minimum can be computed as follows:

Minimize[{ob,c1,c2}, var]

{0, {x → 0, y → 0}}
We can interpret the problem graphically by drawing the plane x + y only in the region: {0 ≤ x ≤ 1, 0 ≤ y ≤ 2}.

Plot3D[x + y, {x, -1, 2}, {y, -1, 2},
RegionFunction -> Function[{x, y}, 0 ≤ x < 1 && 0 ≤ y < 2],
AxesLabel → Automatic, Mesh → None]

The graph clearly shows that the maximum is in (1,2) and the minimum in (0,0). Use the mouse to move the graph around to see it from different perspectives.

In certain situations we’d like the variables to take only integer values. To add this extra requirement in Mathematica, we use the built-in symbol Integers.
Let’s consider the following model in which we want to minimize the objective function (ob) (notice that equality when typing constraints is "=="):

Clear["`Global`*"]
var = {x, y};
ob = x + 2 y;
c1 = - 5 x + y == 7;
c2 = x + y ≥ 26;
c3 = x ≥ 3;
c4 = y ≥ 3;

NMinimize[{ob, c1, c2, c3, c4, var ∊ Integers}, var]

{58., {x → 4, y → 27}}

Since this is a linear problem, we can also use LinearProgramming. As a matter of fact, this command is the most appropriate one for linear problems, especially if the number of variables is large. The syntax is: LinearProgramming [c, m, b]. This function finds the vector x that minimizes the quantity c.x subject to the constraints m.x ≥ b and x ≥ 0. We can limit the values that the variables (some or all of them) can take to just integers with LinearProgramming [..., Integers].

For comparison purposes, let’s solve the same problem using LinearProgramming : x

x+2y→c:{1,2}−5x+y=7→m1:{−5,1}b1:{7,0}x+y≥26→m2:{1,1}b2:{26,1}

Notice that the syntax to indicate the type of constraint is as follows: {b_i, 0} if m_i.x == b_i; {b_i, 1} if m_i.x ≥ b_i and {b_i, -1} if m_i.x ≤ b_i.

LinearProgramming[{1, 2}, {{-5, 1}, {1, 1}}, {{7, 0}, {26, 1}},
{{3, Infinity}, {4, Infinity}}, Integers] // Quiet

{4, 27}

Next we’ll show some examples for nonlinear optimization.

Clear["`Global`*"]
In this problem, both, the objective function that we want to minimize and the constraints are nonlinear.

var = {x, y};
ob = x^2 + (y - 1)^2;
c1 = x^2 + y^2 ≤ 4;
Let’s calculate the global minimum with Minimize:

Minimize[{ob, c1}, var]

{0, {x → 0, y → 1}}
NMinimize can solve this type of problem faster (although in this example we wouldn’t be able notice any difference). We add Chop to remove, when the solution is 0, the spurious values that may appear in the solution when using numerical approximations. Remove Chop from the command below and see what happens.

NMinimize[{ob, c1}, var] // Chop

{0, {x → 0, y → 1.}}
If we are looking for a local minimum we can use FindMinimum. Although it’s not necessary, this function is much faster if we enter an initial point from which to start the search. It’s also a good idea to add Chop.

FindMinimum[{ob, c1}, {{x, 0}, {y, 0}}] // Chop

{0, {x → 0, y → 1.}}
We present the previous problem graphically with Plot3D and ContourPlot, very convenient functions when trying to find optimum points visually. (The option Contours specifies how many contours we want to display; if not specified, Mathematica will choose the optimum number for the given function).

GraphicsRow[{Plot3D[ x^2 + (y - 1)^2, {x, -3, 3},
   {y, -3, 3}, RegionFunction -> Function[{x, y}, x^2 + y^2 ≤ 4],
   AxesLabel → Automatic, Mesh → All], ContourPlot[x^2 + (y - 1)^2,
   {x, -3, 3}, {y, -3, 3}, RegionFunction → Function[{x, y}, x^2 + y^2 ≤ 4],
   Contours → 100, FrameLabel → Automatic]}]

Click inside the output from Plot3D and move the graph around to see the minimum from different angles. In the contour plot, we can see that the minimum corresponds to the coordinates (0, 1). By default, dark colors are associated to small values and light hues are linked to bigger ones. If you place your cursor over the contour limits, you’ll see the value that the function takes at that location. (All the points connected by the same line have the same function value).
Here is another example of nonlinear optimization, in this case for finding the global minimum (Minimize). The objective function is Exp(-x y) with the constraint that x, y ∈ a Circle centered in {0, 0} and with a radius r = 1. We show the result, both numerically and graphically.

m = Minimize[{Exp[-x y], {x, y} ∊ Disk[]}, {x, y}]

{1e,{x→−12,y→−12}}

g = Plot3D[Exp[-x y], {x, y} ∊ Disk[], Axes → True,
Boxed → False, PlotStyle → Opacity@0.6, AxesLabel → Automatic];
Show[g,
Graphics3D[{PointSize@Large, Red, Point[{x, y, m[[1]]} /. m[[2]]]}]]
The following example shows how the optimization algorithm approaches the optimum (minimum) point for f(x,y) in an iterative way.

Clear["Global`*"]

f[x_, y_] := (x - 1)^2 + 2 (y - x^2)^2

Module[{pts}, pts = Last[Last[Reap[Sow[{-1.2, 1, f[-1.2, 1]}];
     FindMinimum[f[x, y],
      {{x, -1.2}, {y, 1}}, StepMonitor :→ Sow[{x, y, f[x, y]}]]]]];
Show[Plot3D[f[x, y], {x, -1.4, 1.4}, {y, -0.6, 2}, Mesh → 5,
   PlotStyle → Automatic, PlotLabel → Style[Row[{"Finding the minimum of ",
       f[x, y], " using iteration"}], 12, Bold], AxesLabel → Automatic],
  Graphics3D[{Orange, PointSize[0.025], Point[Rest[Most[pts]]],
    Thick, Line[pts], Darker@Green, PointSize[0.05],
    Point[Last[pts]], Red, Point[First[pts]]}]]]

11.3.2 Application: Supplying Construction Sites

We’d like to optimize the supply of bricks from warehouses to construction sites. There are two brick warehouses j = {1,2}, where we have in stock 20 and 25 tonnes respectively. From them we supply 6 sites i = {1, 2, 3, 4, 5, 6} with a daily demand in tonnes of: {3, 5, 4, 7, 6, 11}. We assume that the location (latitude and longitude) of the warehouses and sites is:

warehouses = {{7, 1}, {2, 9}};

sites = {{1.25, 1.25}, {8.75, 0.75},
{0.5, 4.75}, {5.75, 5}, {3, 6.5}, {7.25, 7.75}};

We’d like to minimize the total amount of bricks transported per unit of distance, ∑_{i, j} x_i,j d_i,j, meeting the demand without exceeding the stock available in the warehouses.

For that purpose, we have to calculate the minimum of ∑_{i, j} x_i,j d_i,j, where x_i,j corresponds to the quantities transported from the warehouses i:{1, 2} to the sites j:{1, ..., 6} with d_i,j representing the distance between warehouse i and site j.

To calculate the distances between the warehouses and the sites we use the function EuclideanDistance (in the bi-dimensional case, it’s the equivalent to the Pythagorean theorem):

EuclideanDistance[{x1, y1}, {x2, y2}]

√(Abs[x1 - x2]² + Abs[y1 - y2]²
The command below returns triplets, {i, j, d_i,j}, indicating the distance in km between the warehouse i and site j. We use Round[value, 0.1] to round the distance to the nearest first decimal.

Table[{i, j, Round[EuclideanDistance[warehouses[[i]], sites[[j]]], 0.1]},
{i, 1, 2}, {j, 1, 6}]

{{{1, 1, 5.8}, {1, 2, 1.8}, {1, 3, 7.5}, {1, 4, 4.2},
{1, 5, 6.8}, {1, 6, 6.8}}, {{2, 1, 7.8}, {2, 2, 10.7},
{2, 3, 4.5}, {2, 4, 5.5}, {2, 5, 2.7}, {2, 6, 5.4}}}
Therefore, the variables and the objective function, ∑_i,j x_i,j d_i,j, are:

var = {x₁₁, x₁₂, x₁₃, x₁₄, x₁₅, x₁₆, x₂₁, x₂₂, x₂₃, x₂₄, x₂₅, x₂₆};

of = 5.8 x₁₁ + 1.8 x₁₂ + 7.5 x₁₃ + 4.2 x₁₄ + 6.8 x₁₅ +
6.8 x₁₆ + 7.8 x₂₁ + 10.7 x₂₂ + 4.5 x₂₃ + 5.5 x₂₄ + 2.7 x₂₅ + 5.4 x₂₆;
The constraints are: ∑_i x_i,j = demand_j, ∑_j x_i,j = stock_i and x_i,j ≥ 0

constraints = {x₁₁ + x₂₁ == 3 && x₁₂ + x₂₂ == 5 && x₁₃ + x₂₃ == 4 && x₁₄ + x₂₄ == 7 &&
    x₁₅ + x₂₅ == 6 && x₁₆ + x₂₆ == 11 && x₁₁ + x₁₂ + x₁₃ + x₁₄ + x₁₅ + x₁₆ ≤ 20 &&
    x₂₁ + x₂₂ + x₂₃ + x₂₄ + x₂₅ + x₂₆ ≤ 25 && x₁₁ ≥ 0 && x₁₂ ≥ 0 && x₁₃ ≥ 0 && x₁₄ ≥ 0,
  x₁₅ ≥ 0 && x₁₆ ≥ 0 && x₂₁ ≥ 0 && x₂₂ ≥ 0 && x₂₃ ≥ 0 && x₂₄ ≥ 0 && x₂₅ ≥ 0 && x₂₆ ≥ 0};
The solution is:

NMinimize[{of, constraints}, var ]

{149.4, {x₁₁ → 3., x₁₂ → 5., x₁₃ → 0., x₁₄ → 7., x₁₅ → 0.,
x₁₆ → 0., x₂₁ → 0., x₂₂ → 0., x₂₃ → 4., x₂₄ → 0., x₂₅ → 6., x₂₆ → 11.}}

This result indicates that the optimum solution, 149.4 tonnes ×km, is obtained by transporting 3 tonnes from warehouse 1 to site 1, 5 tonnes from warehouse 1 to site 2 and so on.

In “Optimal Transport Scheduling”, a demonstration by Yifan Hu, you can see a dynamic visualization of the resolution of this problem:

http://demonstrations.wolfram.com/OptimalTransportScheduling

If in the demonstration you use 20 and 25 tonnes, you’ll see that the solution is the same (there’s actually a small difference due to rounding).

11.3.3 Application: Portfolio Optimization

A bank commissions a consulting firm to analyze the most profitable way to invest €10,000,000 for two years given the options in the table below. The table includes the expected rate of return (bi-annual) and associated risk.

Product (i)	Return (%), b_i	Risk (%) r_i
Mortgages	9	3
Mutual funds	12	6
Personal loans	15	8
Commercial loans	8	2
Certificates/Bonds	6	1

The capital not invested in any of the products will be placed in government bonds (assumed to be riskless) with a bi-annual rate of return of 3%. The objective of the consulting firm is to allocate the capital to each of the products to meet the following goals:

(a) Maximize the return per € invested.

(b) Keep the possibility of loss to a maximum of 5% of the total amount invested.

(d) Allocate to mutual funds and personal loans an amount no larger than the one invested in mortgages.

Variables: x_i is the percentage of capital invested in product i. The amount placed in government bonds can be considered as a new product with an expected return, in percentage, b₆ = 3 and risk r₆ = 0. Therefore:

var = {x₁(*Mortgages*), x₂(*Mutual funds*),
x₃(*Personal loans*), x₄ (*Commercial loans*), x₅
(*Certificates/Bonds*), x₆(*Government debt*)};
Objective function to maximize: Σi=16bixi.

of = 9 x₁ + 12 x₂ + 15 x₃ + 8 x₄ + 6 x₅ + 3 x₆;
Constraint 1: All the money is invested: Σi=16xi=1.

c1 = x₁ + x₂ + x₃ + x₄ + x₅ + x₆ == 1;
Constraint 2: The average risk is R=Σi=16rixi/Σi=6xi≤5 as r₆ = 0 then:

c2 = (3 x₁ + 6 x₂ + 8 x₃ + 2 x₄ + x₅) / (x₁ + x₂ + x₃ + x₄ + x₅) ≤ 5

(3 x₁ + 6 x₂ + 8 x₃ + 2 x₄ + x₅) / (x₁ + x₂ + x₃ + x₄ + x₅) ≤ 5
Constraint 3: At least 20% of the capital has to be invested in commercial loans (x₄ ≥ 0.2).

c3 = x₄ ≥ 0.2;
Constraint 4: The percentage invested in mutual funds (x₂) and personal loans (x₃) cannot be bigger than the one invested in mortgages (x₁).

c4 = x₂ + x₃ ≤ x₁;
Constraint 5: No percentage invested can be negative (x_i ≥ 0).

c5 = Map[# ≥ 0 &, var];
Using NMaximize, Mathematica tells us that we’d get a return of 11.2% investing 40% in mortgages, 40% in personal loans and 20% in commercial loans.

sol = NMaximize[{of, c1, c2, c3, c4, c5}, var]

{11.2, {x₁ → 0.4, x₂ → 0., x₃ → 0.4, x₄ → 0.2, x₅ → 0., x₆ → 0.}}

"€" <> ToString@AccountingForm[10 000 000 sol[[1]] /100, DigitBlock → 3]

€1,120,000.

11.3.4 Application: Calculating Cod Consumption

Clear["Global`*"]

Many problems in economics involve Cobb–Douglas functions: F(x, y) = A x^a y^b with x, y ≥ 0 and A, a, b constants. Let’s see one example:

The cod consumption in certain region follows the function c(p, i) = 1.3 p^-1.5 i^0.01 where p is the price per kg and i is the daily income in €/day. The previous expression is valid for daily incomes between €20 and €200 per day and the price per kg varies between €4 and €12. What would be the maximum and minimum consumption?

We can proceed as usual defining the variables, objective function and constraints:

c[p_, i_] = 1.3 p^-1.05 i^0.01;

NMinimize[{c[p, i], 4 ≤ p ≤ 12 && 20 ≤ i ≤ 200}, {p, i}]

{0.0985856, {p → 12., i → 20.}}

NMaximize[{c[p, i], 0.4 ≤ p ≤ 1.2 && 20 ≤ i ≤ 200}, {p, i}]

{3.58749, {p → 0.4, i → 200.}}
We can verify the solution visually using Plot2D and ContourPlot. If we place the cursor on the second graph and move it around, we can see the different values of the function c(r, t). This way we can check that the upper left corner corresponds to the largest value (maximum) and the lower right corner to the smallest one (minimum). The boundaries of the curve correspond to the constraints.

GraphicsRow[{Plot3D[c[p, i], {p, 4, 12}, {i, 0, 200},
   RegionFunction -> Function[{p, i}, 4 <= p <= 12 && 20 <= i <= 200],
   AxesLabel -> {"Price", "Income", "c"}, Mesh -> None],
  ContourPlot[c[p, i], {p, 4, 12}, {i, 0, 200}, RegionFunction ->
    Function[{p, i}, 4 <= p <= 12 && 20 <= i <= 200], Mesh -> None,
   Contours -> 100, FrameLabel -> Automatic, PlotLegends → Automatic]},
Frame -> All, AspectRatio → 1/3, ImageSize → {600, 300}]

11.3.5 Application: Transportation Optimization using Linear Programming

The department store chain Carrefa would like to outsource the supply of milk cartons for its centers in Barcelona, Madrid, Seville and Valladolid. For that purpose it receives offers from three suppliers: Pascualo, Simone, and Pulova. The prices per carton are as follows: Pascualo €3.00, Simone €2.80 and Pulova €2.70. The maximum number of cartons that each supplier can send daily is shown in the table below, under the column “Maximum supply”. The daily demand for each center is in the “Demand” row. The rest of the information in the table relates to the transportation costs per carton, in €.

images

Carrefa would like to know how much to order from each supplier to minimize the total cost.

Clear["Global`*"]

This is the type of problem that can be solved with LinearProgramming. For didactic purposes we’re going to compare the notation in NMinimize with the one in LinearProgramming. In practice this is not necessary. As a matter of fact, the easiest way to solve the problem is by using the LinearProgramming syntax.

With NMinimize, the variables can be defined as follows: x_ij = number of cartons from supplier i, i:{1, 2, 3}, with {1 = Pascualo, 2 = Simone, 3 = Pulova}, to center j: {1, 2, 3, 4}, with {1 = Barcelona, 2 = Madrid, 3 = Seville, 4 = Valladolid}. With LinearProgramming (LP), it’s not necessary to define the variables in an explicit way.

{x₁₁, x₁₂, x₁₃, x₁₄, x₂₁, x₂₂, x₂₃, x₂₄, x₃₁, x₃₂, x₃₃, x₃₄} is the same as {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}.

To define the objective function, we have to group terms by supplier. For example: the total cost of sending x units from Pascualo to Barcelona (x₁₁) would be the transportation cost plus the purchasing price: (0.5 + pa). We use the same approach for calculating the rest of the variables:

{pa, sim, pul} = {3, 2.8, 2.7}; LinearProgramming

of = {(0.5 + pa), (0.6 + pa), (0.7 + pa), (0.5 + pa),
(0.5 + sim), (0.6 + sim), (0.8 + sim), (0.5 + sim),
(0.4 + pul), (0.5 + pul), (0.9 + pul), (0.7 + pul)};
The constraint c1, x₁₁ + x₁₂ + x₁₃ + x₁₄ ≤ 2500, representing the maximum number of number of cartons of milk that we can receive from Pascualo is written as (the values of those variables not related to the constraint should be 0):

m1 = {1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0}; b1 = {2500, -1};

where -1 is equivalent to "≤"; 0 to "=="; and 1 to "≥".
The constraint c2, x₂₁ + x₂₂ + x₂₃ + x₂₄ ≤ 1700, related to the maximum capacity of Simone is:

m2 = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0}; b2 = {1700, -1};
The constraint c3, x₃₁ + x₃₂ + x₃₃ + x₃₄ ≤ 1500, dealing with the supply coming from Pulova:

m3 = {0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1}; b3 = {1500, -1};
The constraint c4, x₁₁ + x₂₁ + x₃₁ 800, indicating the demand from Barcelona:

m4 = {1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0}; b4 = {800, 0};
The constraint c5, x₁₂ + x₂₂ + x₃₂ 1100, referring to the demand from Madrid:

m5 = {0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0}; b5 = {1100, 0};
The constraint c6, x₁₃ + x₂₃ + x₃₃ == 500, about Seville:

m6 = {0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0}; b6 = {500, 0};
And finally the constraint c7, x₁₄ + x₂₄ + x₃₄ == 300, to make sure that the needs of Valladolid are being satisfied:

m7 = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1}; b7 = {300, 0};
Now we can solve the problem (LinearProgramming assumes that all the variables are nonnegative so we don’t have to include those constraints).

sol = LinearProgramming[of,
{m1, m2, m3, m4, m5, m6, m7}, {b1, b2, b3, b4, b5, b6, b7}]

{0., 0., 0., 0., 400., 0., 500., 300., 400., 1100., 0., 0.}
LinearProgramming returns the optimal value for each variable. If we’d like to calculate the minimum total transportation cost, we just need to multiply the objective function by the solution (we use "." for matrix multiplication).

of.sol

8870.

11.3.6 Application: Analyzing Enriched Uranium

A key aspect in the operation of nuclear power plants is the cost of the fuel. Virtually all nuclear reactors use uranium as their main source of energy input. The minerals containing it are widely distributed in nature, but profitable deposits under current economic conditions are concentrated in a limited number of countries. According to estimates, these deposits will last for another 70 to 80 years which reduces the amount of money invested in finding new sources. Fast reactors that reuse uranium would extend the availability of the nuclear fuel to thousands of years. However, the cost and the immaturity of the technology don’t justify the development of this type of reactor for the time being. Nevertheless, the current ones, although inefficiently, can reuse part of the spent fuel. Some power plants already do that.

Uranium found in nature mostly consists of three isotopes. The command below displays their respective abundance (fraction of atoms of an isotope in the element). We can see that almost everything is either U238 or U235:

DeleteCases[Transpose[{IsotopeData["U", "Symbol"],
IsotopeData["U", "IsotopeAbundance"] // Chop}], {_, _Integer}]
{{²³⁴U, 0.000054000|, {²³⁵U, 0.007204}, {²³⁸U, 0.992742}}

Most nuclear reactors use uranium enriched in the U–235 isotope. The proportion of this isotope in heavy metals can reach up to 5%. The process to enrich uranium consists of two steps: turning natural uranium into gas (UF6) and then putting the gas into centrifuge tubes (there are other enrichment techniques but this one is the most commonly used). The centrifuge tubes are cylinders that spin at very high speeds pulling the heavier U–238 into the walls of the cylinder and leaving the lighter U–235 closer to the center where it can be sucked out. The result is the extraction of uranium with a higher proportion of U–235 than the one found in nature.

Summarizing: In the enrichment process we start with natural uranium, that has a mass fraction X_f = 0.00711 (equivalent to 0.711% of Uranium–235, or 0.711% in isotopic abundance) and we get enriched uranium with a fraction X_p of U–235 higher than X_f. The process creates depleted uranium as a by-product with a proportion X_w of U–235, lower than the natural one.

The equation below returns the initial number of kilograms needed (kgU of natural uranium, known as feed) to obtain 1 kgU with enrichment X_p, generating as a by-product depleted uranium with a fraction X_w of U–235:

factorfeed=Xp−XwXf−Xw

From this equation we can see that to get 1 kgU with an enrichment X_p, the initial feed depends on the desired content, X_w, of the depleted uranium. The value of X_w expressed as a percentage of U–235 is known as the “tail”. The lower the required X_w, the lower the feed needed. In principle, it may seem that the best outcome would be to obtain a value of X_w as low as technologically possible. However, there’s a problem: the energy consumption needed to generate a target X_w increases nonlinearly as X_w gets smaller. The process is described next.

The method for calculating the energy consumption uses an internationally recognized measurement unit named SWU (Separative Work Unit). This unit is defined as the minimum energy required to obtain 1 kgU enriched by a fraction X_p from uranium with an isotopic fraction X_f (usually natural uranium), generating as a by-product uranium with an isotopic fraction X_w. The actual calculation is done with the following equation:

factorSWU=((2Xp−1)Log[Xp1−Xp]−(2Xw−1)Log[Xw1−Xw])− ((2Xf−1)Log[Xf1−Xf]−(2Xw−1)Log[Xw1−Xw])(Xp−Xw))/(Xf−Xw)

(Log is the natural logarithm, that is, in base e

The previous equation doesn’t refer to the real energy consumption that depends on the technology used. It’s based on the thermodynamical notion of consumption and there’s an agreement to use it regardless of the real consumption.

The two previous equations can be put together in the following function that calculates the feed and SWU necessary to get 1 kgU enriched by 100 X_p and a fixed percentage 100 X_w of depleted U–235. In this function we use W and P to refer to X_w and X_p respectively. F represents X_f and its value is 0.00711 or 0.711% (natural uranium).

swufeed[ P1_, W1_] =
  Module[{F, swu, feed, P, W}, P = P1/100; W = W1/100; F = 0.00711;
swu = ((2*P - 1) *Log[P/ (1 - P)] - (2*W - 1) *Log[W/ (1 - W)]) -
(((2*F - 1) *Log[F/ (1 - F)] - (2*W - 1) *Log[W/ (1 - W)]) * (P - W)) /
      (F - W);
feed = (P - W) / (F - W); Thread[{"factorSWU", "factorFeed"} ->
      {swu, feed}]];
In the next example we calculate the feed and SWU needed to obtain 5% enriched uranium with a tail of 0.3%.

swufeed[ 5, 0.3]

{factorSWU → 7.19832, factorFeed → 11.4355}
For a given tail, in the command below we use 0.2%, the function calculates the SWU required for 1 kgU with enrichment enr:

func [enr_]="factorSWU"/. swufeed[enr,0.2]−6.18776+258.096(−0.002+enr100)+(−1+enr50) Log[enr100(1−enr100)]−6.18776+258.096(−0.002+enr100)+(−1+enr50) Log[enr100(1−enr100)]
We fit the previous function to a line in the range enr 1.0% - 3.5%.

linear1[x_] = Fit[Table[{i, func[i]}, {i, 1.0, 3.5, 0.5}], {1, x}, x]

-1.76531 + 2.02534 x
Now we compare the exact calculation and its numerical approximation. We can see that both lines, although close to each other, are different.

Plot[{linear1[x], func[x]}, {x, 1.0, 3.5},
AxesLabel → {"Enrichment level (%)", "SWU"}]

For the uranium buyer, the setting of the tail is very important and that choice will be determined by the cost of the feed and the SWU.

The enriched uranium cost can be computed using the following expression:

priceEu=factorfeed × pFeed +factorSWU × pSWU, with pFeed being the price of the feed and pSWU the price of the SWU.

The next function calculates the price of 1 kgU for a given tail, enrichment level and feed and SWU prices:

priceEu[enr_, tail_, pFeed_, pSWU_] :=
pFeed "factorFeed" + pSWU "factorSWU" /. swufeed[ enr, tail]

Using the import capabilities discussed in Chapter 2, we can get the spot feed (natural UF6 in kgU) and SWU prices from http://www.uxc.com/review/UxCPrices.aspx. If we open the link, the information will be inside the table: Ux Month-End Spot Prices as of MM, YYYY. The command below will work as long as the website has not been modified since the last time it was accessed (updated on December 31, 2016):

priceuranium = Import[
"http://www.uxc.com/review/uxCPrices.aspx", "Data"][[2, 2, 3 ;; 8]];

Note: To avoid potential problems we have also copied the output directly in Mathematica:

priceuranium =
  {{"U 3 O 8 Price (lb)", "$20.25", "( +2.00)", "€19.35", "( +1.91)"},
   {"NA Conv. (kgU)", "$5.85", "( Unch.)", "€5.59", "( Unch.)"},
   {"EU Conv. (kgU)", "$6.40", "( Unch.)", "€6.12", "( Unch.)"},
   {"NA UF 6 Price (kgU)", "$58.00", "( +4.60)", "€55.43", "( +4.40)"},
   {"NA UF 6 Value § (kgU)", "$58.76", "( +5.23)",
    "€56.15", "( +5.00)"}, {"EU UF 6 Value § (kgU)",
    "$59.31", "( +5.23)", "€56.68", "( +5.00)"}};

We can even choose a specific value for later use. The next command extracts the price of the feed (EU UF6 value). Nevertheless, we copy it so that the reader can reproduce the example using the same data.

{priceuranium[[5, 2]], priceuranium[[6, 2]]}

{$58.76, $59.31}

{pFeed, pSWU} = {58.76, 59.31};
Given the spot prices for the feed and SWU, the cost of 1 kgU enriched by 5% and with a 0.3% tail will be:

priceEu[5, 0.3, pFeed, pSWU]

1098.88
The graph below displays the price as a function of the tail. We can clearly see that there’s a minimum economic value associated to the tail.

ListPlot[Table[{i, priceEu[5, i, pFeed, pSWU]}, {i, 0.10, 0.40, 0.025}],
AxesOrigin → {0, 2500}, Joined → True, AxesLabel → {"Tail(%)", "Price"}]
The following function calculates the optimum purchasing prices for different enrichment levels. Notice that for the given feed and SWU prices, the optimum tail doesn’t depend on the enrichment target.

Table[{enr, NMinimize[{priceEu[enr, tail, pFeed, pSWU], 0.15 . tail . 0.5},
tail]}, {enr, 1, 5}]

{{1, {114.273, {tail → 0.22767}}},
{2, {335.907, {tail → 0.22767}}}, {3, {575.243, {tail → 0.22767}}},
{4, {821.986, {tail → 0.22767}}}, {5, {1072.88, {tail → 0.22767}}}}

In practice, the tail value is rounded to two decimals.

The next example shows how to calculate the cost of the enriched uranium needed to generate 1 kWh in European Union.

The normal annual consumption of a 1,000 MW power plant, the most common one, is of 20,000 kgU 4.5% enriched. Using the optimum tail from the previous data, the cost of this uranium expressed in euros would be:

eurodollar = FinancialData["EUR/USD"]

1.0601

annualcomsumption = eurodollar 20 000 priceEu[4.5, 0.23, pFeed, pSWU] "euros"

2.00794×10⁷ euros
A 1,000 MW power plant working non-stop for one year will generate the following amount of energy in kWh (kilowatt-hour):

annualkwh = 1 000 000 × 24 × 365 "kWh"

8 760 000 000 kWh
Therefore, the cost per kWh would be:

annualcomsumption /annualkwh

0.00229217 euroskwh
The previous cost refers exclusively to the enriched uranium. To that number we must add the cost of its production, transportation and insurance among others. In total, around 20% of the uranium price. Therefore, the cost of fuel per kWh expressed in euro cents would be:

nuclearfuelkwh = 100 × 1.2 annualcomsumption/annualkwh/"euros" "cents"

0.27506 centskwh

This figure doesn’t reflect the true cost of nuclear kWh, that is influenced by additional factors, especially the initial construction investment.

The cost of the electricity in Europe is higher than 10 c/kWh. If we compare it with the nuclear fuel cost, we can see that the latter represents only a small portion of the former. Therefore, big swings in fuel costs will have a small impact in the cost per nuclear kWh. In the case of other fuels such as natural gas, the most commonly used one for generating electricity, their costs can represent up to 80% of the total production expense.

11.4 The Shortest Path Problem

11.4.1 The Traveling Salesman Problem

There are problems in mathematics that look very specialized and easy to solve at first sight but in reality they are very complicated and their resolution can have applications in a wide variety of fields. One of the most famous examples is the Traveling Salesman Problem (TSP). This problem can be defined as follows: A traveler would like to visit N cities starting and finishing in the same arbitrarily chosen one, minimizing the total distance traveled and passing through each city only once. If you think about it, there are many real-world problems that are actually equivalent to the TSP: route optimization (transportation, logistics), optimum network layout (trains, roads, electricity), even the connections inside microprocessors to minimize calculation time.

A graphical demonstration of the problem created by Jon McLoone can be found in: (http://demonstrations.wolfram.com/TravelingSalesmanProblem/)

Mathematically, the problem consists of finding a permutation P = {c₀, c₁,..., c_n-1} such that dP=Σi=1n−1d[ci,ci+1] would be a minimum, with d_ij representing the distance from city i to city j.

The seemingly simplest solution would be to enumerate all the possible routes, calculate the distance traveled for each route and select the shortest one. The problem with this approach is that computing time increases dramatically when the number of cities gets bigger. For example, a computer would take 3 seconds to solve the problem if there were 10 cities, a bit more than half a minute for 11 cities and almost 80,000 years! to solve the problem for 20 cities.

Conclusion: This solving method is not feasible. Therefore, several alternative algorithms have been developed that use a variety of strategies to simplify the problem (Something similar to computer chess programs that instead of calculating all possible moves, use different criteria to choose the most appropriate ones).

For this type of problem, Mathematica has the function:

FindShortestTour[{e1, e2, ...}]. This function takes a list of coordinates: {e1, e2, ...} and attempts to find the optimum ordering so that the total distance traveled is minimized and all the points have been visited just once. The output returns both, the total distance and the order of traveling.

Example: Given a set of points p = {1, 2, 3, 4, 5, 6} each with coordinates{x,y}:{{4, 3},{1, 1}, {2, 3},{3, -5},{-1, 2},{3, 4}}, find the shortest path that will visit all the points only once.

d = {{4, 3}, {1, 1}, {2, 3}, {3, -5}, {-1, 2}, {3, 4}};

{dist, order} = FindShortestTour[d]

{22+5+310+65,{1,4,2,5,3,6,1}}

We sort the points following the order from the output above and represent them graphically:

d[[order]]

{{4, 3}, {3, -5}, {1, 1}, {-1, 2}, {2, 3}, {3, 4}, {4, 3}}

Graphics[Line[%]]

11.4.2 A Tour Around South American Cities

The next example is about organizing a trip to several cities in South America. In what order should we visit them to minimize the total distance traveled?

The cities we’re visiting are the following (we indicate both, the cities and their respective regions and countries to avoid ambiguity):

cities={Entity["City", {"Asuncion", "Asuncion", "Paraguay"}],
Entity["City", {"Bogota", "DistritoCapital", "Colombia"}],
Entity["City", {"RioDeJaneiro", "RioDeJaneiro", "Brazil"}],
Entity["City", {"BuenosAires", "BuenosAires", "Argentina"}],
Entity["City", {"Caracas", "DistritoCapital", "Venezuela"}],
Entity["City", {"LaPaz", "LaPaz", "Bolivia"}],
Entity["City", {"Lima", "Lima", "Peru"}],
Entity["City", {"Montevideo", "Montevideo", "Uruguay"}],
Entity["City", {"Quito", "Pichincha", "Ecuador"}],
Entity["City", {"Santiago", "Metropolitana", "Chile"}]};
Now we can calculate the best visiting order:

order = Last[FindShortestTour[GeoPosition[cities]]]

{1, 8, 4, 10, 6, 7, 9, 2, 5, 3, 1}
Finally, we can display the route on a map:

GeoListPlot[cities[[order]], Joined → True]

11.4.3 Moon Tourism

It’s the year 2035 and we are planning to organize sightseeing tours around the Moon. We want to find the shortest route for visiting all the Apollo landing sites.

“Apollo landings” is one of the available entity classes that have been included in the function MannedSpaceMissionData.

MannedSpaceMissionData["Classes"]
The names of the manned landing missions are:
Notice that Apollo 13 is missing. Let’s find out why:

{intended to be the third manned mission to land on the moon,
mid-mission onboard explosion forced landing to be aborted,
returned safely to Earth}
Create a map showing the locations of the landings:

GeoListPlot[missionnames, GeoRange → All,
GeoProjection → "Orthographic", GeoLabels → True, LabelStyle → White]
To extract the coordinates of those landings:

landingcoords = Values[apollolandings] // Flatten
Next, we find out that the shortest route starts where Apollo 11 landed:

order = Last[FindShortestTour[landingcoords]]

{1, 6, 4, 2, 3, 5, 1}
Finally, we create a plot showing the path of the tour:

GeoListPlot[landingcoords[[order]], Joined → True, GeoLabels → Automatic]

11.5 Optimum Flows

In this example we want to maximize the total number of freight cars that can be sent from Vancouver to Winnipeg given the following network characteristics:

The code generating graph g is not shown.

images

The previous illustration has been created using the Graph function (for further details see Section 4.4). For simplifying purposes, let’s just focus on 3 cities: { (1) Vancouver, (2) Calgary, (3) Edmonton}. The flow directions (i→j) and capacities (EdgeCapacity) are: {From (1) Vancouver to (2) Calgary, 13, from (1) Vancouver to (3) Edmonton, 4, from (3) Edmonton to (2) Calgary, 10 and from (2) Edmonton to (3) Calgary, 4}. We display that information using DirectedEdge and with the help of VertexLabels and Placed position the labels in the desired locations. The style is defined using Style [#, “Style”]&. You can use the documentation to learn more about the syntax for each function.

g1 = Graph[{
    "Vancouver", "Calgary", "Edmonton"}, {{{
     1, 2}, {1, 3}, {2, 3}, {3, 2}}, Null}, {
    EdgeCapacity -> {13, 16, 4, 10},
     EdgeLabels -> {
       DirectedEdge["Vancouver", "Calgary"] -> Placed[13, {0.5, {2, 1}}],
      DirectedEdge["Vancouver", "Edmonton"] -> Placed[16, {0.5, {-1.5, 2}}],
    DirectedEdge["Edmonton", "Calgary"] -> Placed[10, {0.5, {-0.2, 1.5}}],
     DirectedEdge["Calgary", "Edmonton"] -> Placed[4, {0.5, {-0.7, 1.5}}] },
    GraphStyle -> "BasicBlack", VertexLabels -> {
     "Edmonton" -> Placed["(3)Ed.", Above, Style[#, Red] & ],
     "Vancouver" -> Placed["(1)Vanc.", After, Style[#, Blue] & ],

     "Calgary" -> Placed["(2)Cal.", Above, Style[#, FontFamily -> "Helvetica"] & ]},
    VertexSize -> {Small}}, ImagePadding → All, ImageSize → Small]
In the next example we add VertexCoordinates to specify the coordinates for the center of the vertices. We also modify several options to make it easier to overlay the graph on a map of Canada later.

Graph[{
    "Vancouver", "Calgary", "Edmonton"}, {{{
     1, 2}, {1, 3}, {2, 3}, {3, 2}}, Null}, {
    EdgeCapacity -> {13, 16, 4, 10},
     EdgeLabels -> {
      DirectedEdge["Edmonton", "Calgary"] -> Placed[10, {0.5, {-0.2, 1}}],
      DirectedEdge["Vancouver", "Edmonton"] -> Placed[16, {0.5, {1, 0}}],
     DirectedEdge["Vancouver", "Calgary"] -> Placed[13, {0.6, {1, 1.5}}]},
     GraphStyle -> "BasicBlack", ImagePadding -> {{0, 40}, {0, 10}},
     VertexCoordinates -> {{-0.253, 0.112}, {-0.1419, 0.1239}, {-0.122,
     0.1735}}, VertexLabels -> {
     "Edmonton" ->
     Placed["Edmonton", Above, Style[#, FontFamily -> "Helvetica"] & ],
     "Vancouver" ->
     Placed["Vancouver", {{1.5, -1}, {0, 0}},
      Style[#, FontFamily -> "Helvetica"] & ],
     "Calgary" ->
     Placed["Calgary", Below, Style[#, FontFamily -> "Helvetica"] & ]},
    VertexSize -> {Small}}]
To find the optimum flow we use the function FindMaximumFlow.

F = FindMaximumFlow[g, "Vancouver", "Winnipeg", "OptimumFlowData"]

F["FlowValue"]

23
We can see the optimum flow per edge as follows:

Text@Grid[Prepend[{#, F[#]} & /@ F["EdgeList"],
{Style["Railroad", Bold], Style[ "Edge Flow", Bold]}]]
Here we display the result on top of the map of Canada, downloaded using CountryData (for further details see Chapter 5).

Panel[Show[{CountryData["Canada", "Shape"], F["FlowGraph"]},
PlotRange → {{-0.31, 0.083}, {0.04, 0.25}}]]

11.6 Additional Resources

The Wolfram Finance Platform: http://www.wolfram.com/finance-platform

Mathematica's finance related functions:

http://reference.wolfram.com/mathematica/guide/Finance.html

Financial engineering: http://www.wolfram.com/solutions/industry/financial-engineering-and-mathematics

Economics: http://www.wolfram.com/solutions/industry/economics/

Financial risk management: http://www.wolfram.com/solutions/industry/financial-risk-management/

Constrained optimization tutorial: tutorial/ConstrainedOptimizationOverview