Serving the Quantitative Finance Community

 
User avatar
Cuchulainn
Posts: 20254
Joined: July 16th, 2004, 7:38 am
Location: 20, 000

Re: exp(5) = [$]e^5[$]

April 18th, 2017, 10:48 am

Instead of Newton, bisection, analytical approximations (too many assumptions that can break during stormy weather) we do least-squares minimisation of the unimodal(!) function [$](log(x) - 5.0)^2[$]  using the (most?) efficient three-point interval search on [32,243] with the least number of function evaluations. Does it do better than other algos for iv?

Here is the output (tol = 1.0e-5). Nice thing is the answer always lies in an interval.

Minimum: x, f(x) [84.75,190.25], 0.187801
Minimum: x, f(x) [111.125,163.875], 0.0467709
Minimum: x, f(x) [137.5,163.875], 0.00782745
Minimum: x, f(x) [144.094,157.281], 0.00212025
Minimum: x, f(x) [144.094,150.688], 0.000551828
Minimum: x, f(x) [147.391,150.688], 0.000139543
Minimum: x, f(x) [147.391,149.039], 3.27546e-05
Minimum: x, f(x) [147.803,148.627], 9.52942e-06
Minimum: x, f(x) [148.215,148.627], 1.93004e-06
Minimum: x, f(x) [148.318,148.524], 4.84548e-07
Minimum: x, f(x) [148.369,148.472], 1.23176e-07
Minimum: x, f(x) [148.395,148.447], 3.28335e-08
Minimum: x, f(x) [148.395,148.421], 8.72959e-09
Minimum: x, f(x) [148.408,148.421], 1.95913e-09
Minimum: x, f(x) [148.411,148.418], 5.47337e-10
Minimum: x, f(x) [148.411,148.414], 1.22007e-10
Minimum: x, f(x) [148.412,148.414], 3.37671e-11
Minimum: x, f(x) [148.413,148.414], 7.74401e-12
Minimum: x, f(x) [148.413,148.413], 2.22907e-12
Minimum: x, f(x) [148.413,148.413], 4.62369e-13
Minimum: x, f(x) [148.413,148.413], 1.17684e-13
Minimum: x, f(x) [148.413,148.413], 3.15131e-14
Minimum: x, f(x) [148.413,148.413], 8.20017e-15
Minimum: x, f(x) [148.413,148.413], 1.90432e-15
Bracket  [148.413,148.413]
Bracket [2.78939e-15,1.01925e-15]

For OP, tol = 1.0e-1 gives

Minimum: x, f(x) [84.75,190.25], 0.187801
Minimum: x, f(x) [111.125,163.875], 0.0467709
Minimum: x, f(x) [137.5,163.875], 0.00782745
Minimum: x, f(x) [144.094,157.281], 0.00212025
Minimum: x, f(x) [144.094,150.688], 0.000551828
Minimum: x, f(x) [147.391,150.688], 0.000139543
Minimum: x, f(x) [147.391,149.039], 3.27546e-05
Minimum: x, f(x) [147.803,148.627], 9.52942e-06
Minimum: x, f(x) [148.215,148.627], 1.93004e-06
Minimum: x, f(x) [148.318,148.524], 4.84548e-07
Minimum: x, f(x) [148.369,148.472], 1.23176e-07
Bracket  [148.369,148.472]
Bracket [8.70207e-08,1.59332e-07]
 
User avatar
Cuchulainn
Posts: 20254
Joined: July 16th, 2004, 7:38 am
Location: 20, 000

Re: exp(5) = [$]e^5[$]

April 21st, 2017, 9:46 am

We all know by now that the answer is 148.413 and that the function to be minimised is simple. For more complicated functions things can become more tricky and then we want to make C++ look more like mathematics (functional programming languages do it before breakfast). C++11 helps a bit by lambda functions so that you can build a vector space of functions to assemble functions from simpler ones. You would kind of expect that a library builder out there has done this (knowledge of maths is needed).
// TestHOF101.cpp
//
// Testing lambda functions; build complex functions using algebra. 
// A vector space of functions.
// In  this case compute x - e^5 = 0 using least squares and Brent's method
// x == 148.413
//
// DD

#include <functional>
#include <cmath>
#include <iostream>
#include <boost/math/tools/minima.hpp>

template <typename T>
 std::function<T (T)> operator * (std::function<T (T)>& f, std::function<T (T)>& g)
{ // Addition

 return [=](T x)
 {
 return  f(x) * g(x);
 };
}

template <typename T>
 std::function<T (T)> operator - (std::function<T (T)>& f, T a)
{ // Addition

 return [=](T x)
 {
 return  f(x) - a;
 };
}



template <typename T>
 std::function<T (T)> log(std::function<T (T)>& f)
{ // log

 return [=](T x)
 {
 return  std::log(f(x));
 };
}
template <typename T>
 std::function<T(T)> exp(std::function<T(T)>& f)
{ // log

 return [=](T x)
 {
 return  std::exp(f(x));
 };
}

double funcGlobal(double x)
{
 return (std::log(x) - 5.0)* (std::log(x) - 5.0);
}

int main()
{

 // Define interval [min, max] where search will take place
 double min = 5.0;
 double max = 300.0;

 int bits = 50;
 boost::uintmax_t maxIter = 10000;

 // Standard function call
 std::pair<double, double> result = boost::math::tools::brent_find_minima(funcGlobal, min, max, bits, maxIter);
 std::cout << "Abscissa, value f(x) : " << result.first << ", " << result.second << std::endl;

 // Now use lambda and algebra; assemble functions incrementally
 std::function<double(double)> f1 = [](double x) { return std::log(x); };
 decltype(f1) f2 = f1 - 5.0;
 decltype(f2) funcHOF = f2*f2;

 // Stress testing
 std::function<double(double)> f0 = [](double x) { return x; };
 decltype(f0) funcExtreme = exp(log(exp(log(f0))));
 std::cout << "? " << funcExtreme(5) << '\n';
 decltype(f0) f3 = log(funcExtreme);
 decltype(f0) f4 = exp(funcExtreme);
 std::cout << "? " << f3(5) << '\n';
 std::cout << "? " << f4(5) << '\n';
 //

 std::pair<double, double> result2 = boost::math::tools::brent_find_minima(funcHOF, min, max, bits, maxIter);
 std::cout << "Abscissa, value f(x) HOF: " << result2.first << ", " << result2.second << std::endl;
 
}
//
 
User avatar
Cuchulainn
Posts: 20254
Joined: July 16th, 2004, 7:38 am
Location: 20, 000

Re: exp(5) = [$]e^5[$]

April 22nd, 2017, 1:21 pm

Until now, we used double precision data types for the minimisation using interval searching. Will the algorithm work with multiprecision types? It would seem yes at the moment. We use 100-digit type and we get 148.413..... to 100 digits that we checked against Boost and that big online calculator.

And the run-time performance doesn't seem to be all that bad. No test done.

The output is

3-point interval Search method
148.4131 5910 2576 6034 2111 5580 0405 5227 9623 4876 6759 3878 9890 4675 2845 1109 1206 4820 9585 7607 9688 4094 5989 90211
 
Exact (yuge precision)
148.4131 5910 2576 6034 2111 5580 0405 5227 9623 4876 6759 3878 9890 4675 2845 1109 1206 4820 9585 7607 9688 4094 5989 9021 1
 
\-> etc.412 9280 8270 6663 2605 2909 9262 3769 4548 1827 8044 9112 4354 2477 7223 6978 7450 8409 7439 3799 2123 5516 6686 3241 8855 9697 6672 2232 1326 5623 4323 8300 9773 3581 8269
 
 
User avatar
Cuchulainn
Posts: 20254
Joined: July 16th, 2004, 7:38 am
Location: 20, 000

Re: exp(5) = [$]e^5[$]

September 7th, 2017, 8:16 pm

What if the entire notion that a transition probability even exists is false?
If you'd read PWOQF2 then you'd know that people might quote transition matrices for a period of a year for which there is no infinitesimal-period matrix, because you can't invert e^M (and by "invert" I don't mean one "one over"!). The implication being either the transition matrix being quoted is rubbish or that whoever generates the matrix is making some sophisticated assumption about time dependence. (No, I don't believe the latter either!) Yet another reason why the proper defintion of exp(M) is important and the other one is pretty pointless!
P
One way to persuade infidels is to define [$]e^A[$] as the solution of a matrix ODE.

[$]dB/dt = AB, B(0) = I[$]

The solution is [$]B(t) = e^{At}[$]. So plug in  [$] t = 1[$] and Bob's your uncle.  

Now, the next step is to approximate the matrix ODE by a numeric ODE solver, easy in C++.
 
User avatar
Cuchulainn
Posts: 20254
Joined: July 16th, 2004, 7:38 am
Location: 20, 000

Re: exp(5) = [$]e^5[$]

September 7th, 2017, 8:21 pm

Here's the C++ code in Boost odeint. Few lines of code.


using value_type = double;
typedef boost::numeric::ublas::matrix<value_type> state_type;

class MatrixOde
{
private:
 // dB/dt = A*B, B(0) = C;
 boost::numeric::ublas::matrix<value_type> A_;
 boost::numeric::ublas::matrix<value_type> C_;
public:
 
 MatrixOde(const boost::numeric::ublas::matrix<value_type>& A, const boost::numeric::ublas::matrix<value_type>& IC)
 : A_(A), C_(IC) {}

    void operator()( const state_type &x , state_type &dxdt , double t ) const
    {
 
        for( std::size_t i=0 ; i < x.size1();++i )
        {
            for( std::size_t j=0 ; j < x.size2(); ++j )
            {
 dxdt(i, j) = 0.0;
 for (std::size_t k = 0; k < x.size2(); ++k)
 {
 dxdt(i, j) += A_(i,k)*x(k,j);
 }
 }
 }
  }

};
 
JGJPR
Posts: 12
Joined: September 7th, 2017, 3:33 pm

Re: exp(5) = [$]e^5[$]

September 10th, 2017, 6:48 pm

Is a quant a calculation engine?
Very strange question to me...
 
User avatar
Cuchulainn
Posts: 20254
Joined: July 16th, 2004, 7:38 am
Location: 20, 000

Re: exp(5) = [$]e^5[$]

September 10th, 2017, 8:29 pm

Is a quant a calculation engine?
Very strange question to me...
Good point. Some want closed-form solution.
 
User avatar
Cuchulainn
Posts: 20254
Joined: July 16th, 2004, 7:38 am
Location: 20, 000

Re: exp(5) = [$]e^5[$]

September 18th, 2017, 8:11 am

Here is another solution; all you need is to use your office wall as a dart board.

[$]\int_{0}^{5} \exp(x) dx = \exp(5) -1 [$]

Then use Monte Carlo Hit or Miss Integration. Draw the curve and throw darts at the wall. How many hits? Maybe get more people to join in? Just be careful when throwing those darts!
 
User avatar
Cuchulainn
Posts: 20254
Joined: July 16th, 2004, 7:38 am
Location: 20, 000

Re: exp(5) = [$]e^5[$]

October 27th, 2017, 8:43 am

The only solution left is tie solve this problem using (a hand-crafted?) ML algorithm?

Can I just look at the diagram and give the answer? 

Image
 
User avatar
Traden4Alpha
Posts: 3300
Joined: September 20th, 2002, 8:30 pm

Re: exp(5) = [$]e^5[$]

October 27th, 2017, 1:15 pm

"Paper and pencil" would let us construct our own slide rule either for multiplication or for direct exponentiation. But it seems like a labor-intensive solution in that: 1) it requires assigning accurate multi-digit decimal values to specific marks on the slide rule; 2) it seems to imply that we must solve the more general problem of [$]e^x[$]. It also might take a lot of paper if that if the marks on the paper are accurate to only O(1 mm), then we might need a piece of paper O(14841 mm) long. We can use a roll of adding machine paper to do this.
 
User avatar
Cuchulainn
Posts: 20254
Joined: July 16th, 2004, 7:38 am
Location: 20, 000

Re: exp(5) = [$]e^5[$]

October 29th, 2017, 11:22 am

Here's the C++ code in Boost odeint. Few lines of code.


using value_type = double;
typedef boost::numeric::ublas::matrix<value_type> state_type;

class MatrixOde
{
private:
 // dB/dt = A*B, B(0) = C;
 boost::numeric::ublas::matrix<value_type> A_;
 boost::numeric::ublas::matrix<value_type> C_;
public:
 
 MatrixOde(const boost::numeric::ublas::matrix<value_type>& A, const boost::numeric::ublas::matrix<value_type>& IC)
 : A_(A), C_(IC) {}

    void operator()( const state_type &x , state_type &dxdt , double t ) const
    {
 
        for( std::size_t i=0 ; i < x.size1();++i )
        {
            for( std::size_t j=0 ; j < x.size2(); ++j )
            {
 dxdt(i, j) = 0.0;
 for (std::size_t k = 0; k < x.size2(); ++k)
 {
 dxdt(i, j) += A_(i,k)*x(k,j);
 }
 }
 }
  }

};
A sanity check is to take the scalar case du/dt = u. I think we all agree that a scalar can be seen as a 1x1 matrix. So, nothing shocking really but it is OP. 
148.413.
// Sanity check; 1X1 matrix case dB/t = B on (0,5), B(0) = I
	{ // dB/dt = A*B, B(0) = C

		std::size_t N = 1;
		ublas::identity_matrix<value_type> C(N);
		ublas::matrix<value_type> A(N, N);
		A(0, 0) = 1.0;
		ublas::matrix<value_type> B(N, N);

		// Initialise the solution matrix B
		B = C;

		MatrixOde ode(A, C);
		double L = 0.0;
		double T = 5.0;
		double dt = 0.001;

		std::size_t steps = Bode::integrate(ode, B, L, T, dt, write_out);

		std::cout << "Value at time: " << T << '\n';
		std::cout << B << '\n';
	}

 
User avatar
Traden4Alpha
Posts: 3300
Joined: September 20th, 2002, 8:30 pm

Re: exp(5) = [$]e^5[$]

November 2nd, 2017, 4:12 pm

Put $1 in a bank account yielding exactly 2% continuously compounded and then check the account balance at noon in 252 years and 179 days.
 
User avatar
Cuchulainn
Posts: 20254
Joined: July 16th, 2004, 7:38 am
Location: 20, 000

Re: exp(5) = [$]e^5[$]

November 5th, 2017, 11:21 am

One glaring omission is that no one has mentioned how [$]e[$] came to be. Everyone blames Euler (some blame APL) but its roots are probably due to a 17th century banker or trader who stumbled on compound interest:

[$]\displaystyle\lim_{n\to\infty} (1 + r/n)^{nt} \rightarrow e^{rt}[$]

Then the special case gives us

[$]\displaystyle\lim_{n\to\infty} (1 + 1/n)^{n} == e[$]

There you have it. Even Einstein agrees on [$]e^5[$].



Image
Last edited by Cuchulainn on November 5th, 2017, 11:37 am, edited 2 times in total.
 
User avatar
Cuchulainn
Posts: 20254
Joined: July 16th, 2004, 7:38 am
Location: 20, 000

Re: exp(5) = [$]e^5[$]

November 5th, 2017, 11:30 am

Put $1 in a bank account yielding exactly 2% continuously compounded and then check the account balance at noon in 252 years and 179 days.
Very good answer. I had not read this until this morning. But funnily the same idea crossed my mind yesterday.

So, the question was: How long does it take to earn $148.413 by investing $1 at a fixed rate of 2%?
 
User avatar
Traden4Alpha
Posts: 3300
Joined: September 20th, 2002, 8:30 pm

Re: exp(5) = [$]e^5[$]

November 5th, 2017, 1:52 pm

Put $1 in a bank account yielding exactly 2% continuously compounded and then check the account balance at noon in 252 years and 179 days.
Very good answer. I had not read this until this morning. But funnily the same idea crossed my mind yesterday.

So, the question was: How long does it take to earn $148.413 by investing $1 at a fixed rate of 2%?
Yep!

I also pondered if other exponential processes might be adapted to answer this question. But base-2 doubling requires an inconvenient 5/ln(2) = 7.213475 doublings to get the answer. And looking for the right seed for a Fibonacci process seemed to yield an answer of starting with 7457 and 12067 which after 10 cycles reaches 10000*e^5 with requisite accuracy but it feel like a hack!