SERVING THE QUANTITATIVE FINANCE COMMUNITY

 
User avatar
Cuchulainn
Posts: 51410
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

Re: System of nonlinear equations

June 19th, 2017, 3:01 pm

It could happen to a bishop! Seriously, people can get stuck in a rut and then always have several solutions.

https://en.wikipedia.org/wiki/How_to_Solve_It

I am surprised DE does not converge in your case. DE is very robust. So, it might be bad input or a bug in Matlab , or other?
http://www.datasimfinancial.com

I tell the kids, somebody’s gotta win, somebody’s gotta lose. Just don’t fight about it. Just try to get better.
 
User avatar
Alan
Posts: 8984
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

Re: System of nonlinear equations

June 19th, 2017, 3:17 pm

@Gencho, glad to hear it worked out. (that google link doesn't work at all).

@Daniel, if you still want the paraphrase: the basic model (Hull example), models the equity as a call option -- with the asset value [$]V_t[$] playing the role of the underlying and a call strike equal to the debt value L. In that model, the company goes bankrupt if and only if the asset value is less than the debt at expiration T. But paths where the intermediate asset values (prior to expiration) fall below the strike (=debt) are fine (do not result in bankruptcy) as long as [$]V_T > L[$]. 

The second model tightens up the bankruptcy criterion so that the company is bankrupt if the asset value i) falls below the barrier H at *any* time prior to T  or ii) is below L at expiration. The natural choice is H = L, but you could argue that H < L could make sense, too. In any event, If H is much smaller than L (H << L), the likelihood of the particle reaching the barrier is small, so it's as if the barrier does not exist, and so you should recover the basic model results.  
 
User avatar
Cuchulainn
Posts: 51410
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

Re: System of nonlinear equations

June 20th, 2017, 11:45 am

Gencho wrote:
P.S. I am very surprised, to say at least, that I could not guess myself the following:
1) I could not understand that scaling is so important.
2) I could not even think that the root solver will work worse than constrained optimization. I thought it is the most robust approach.
So again, you really helped me!

Once a given approach becomes accepted then people just tend to use and if necessary tweak it until it works (or not). In the current case the usual suspect is a variation on Newton Raphson.
Of course the "dropping penny" is the realisation [$] f(x) = 0 \Leftrightarrow min f(x)^2.[$]

There was a discussion on this very topic with many (ingenious) solutions. AFAIR minimisaton was the most robust.

viewtopic.php?f=34&t=97812

//
Another robust method to resolve NR is homotopy/continuation/Davidenko (from topology) is to embed f(x) = 0 in a lager space to get H(x,t) = where 0 <= t <= 1 and solve the corresponding ODE.
http://www.lehigh.edu/~wes1/apci/24feb00-s.pdf

Here's code in 1d for Kepler's equation; In 2d you would us the Jacobian. AFAIR Matlab has a ODE solver ode45??
#include <iostream>
#include <vector>
#include <cmath>
#include <functional>

#include <boost/numeric/odeint.hpp>
#include <boost/numeric/odeint/stepper/bulirsch_stoer.hpp>

// The type of container used to hold the state vector 
typedef std::vector<double> state_type;
typedef std::function<double (double)> FunctionType;

namespace Bode = boost::numeric::odeint;
using namespace boost::numeric::odeint;


double p = 2.0;
double c = 10.0;

const double eps = 13.5;
const double M = 0.8;


double Kepler(double y)
{
 double d = y - eps*std::sin(y) - M;
 return d;
}

// NR function stuff
double f(double x)
{

 return Kepler(x);
}
double fd(double x)
{
 
 // Kepler
 return 1.0 - eps*std::cos(x);
}

//double x0 = 138.3377; // far away
//double x0 = 12.3; // far away
double x0 = 2.3000001; // try -20.0 as well

class WriteOutput
{
public:
 std::vector<double> tValues;
 std::vector<double> yValues;

public:
    WriteOutput() 
 : tValues(std::vector<double>()), yValues(std::vector<double>())
 { }
 
 void operator ()( const state_type &x , const double t )
 {
 tValues.push_back(t); yValues.push_back(x[0]); 

 std::cout << "Time and value (FO): " << t << " " << x[0] << std::endl;
 }
 };
// The rhs of x' = f(x) defined as a class 
class Davidenko
{
public:
    Davidenko() { }

    void operator() (const state_type& x, state_type& dHdt, const double t)
    {
  dHdt[0] = -f(x0)/fd(x[0]);
    }
};


void write_out( const state_type& x , const double t )
{
 if (t == 1.0)
 {
 std::cout << "Time and value: " << t << " " << x[0] << std::endl;
 }
}




int main()
{
    namespace Bode = boost::numeric::odeint;

    // Initial condition
    state_type x(1);
 x[0] = x0;

 double L = -1.0;
 double U = 2.0;
 double dt = 1.0e-3;

 // Integration_class, V2 based on Runge Kutta54 Cash Karp stepper (5th order)
 Davidenko ode;
    std::size_t steps = Bode::integrate(ode, x, L, U, dt, write_out);

 std::cout << "Number of steps " << steps << std::endl;
 std::cout << x[0] << ", " << f(x[0])<< std::endl;
 

 WriteOutput wo;
// Choice B: Fast adaptive
 typedef runge_kutta_cash_karp54< state_type > error_stepper_type;
 typedef Bode::controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
 Bode::runge_kutta_fehlberg78< state_type, double> myStepper;
// steps = Bode::integrate_const(myStepper, ode, x, L, U, dt, boost::ref(wo));
// double abs_err = 1.0e-3 , rel_err = 1.0e-3, a_x = 1.0 , a_dxdt = 1.0; double dtODE = 1.0e-3; 

 //   controlled_stepper_type controlled_stepper(default_error_checker< double >( abs_err , rel_err , a_x , a_dxdt ) );
// steps = Bode::integrate_adaptive( controlled_stepper,ode, x , L, U, dtODE, boost::ref(wo) ); 
 
// std::cout << "Steps " << steps << ", " << f(x[0]) << std::endl;
 
 return 0;
}

Last edited by Cuchulainn on June 20th, 2017, 12:25 pm
http://www.datasimfinancial.com

I tell the kids, somebody’s gotta win, somebody’s gotta lose. Just don’t fight about it. Just try to get better.
 
User avatar
Cuchulainn
Posts: 51410
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

Re: System of nonlinear equations

June 20th, 2017, 11:59 am

Alan wrote:
@Gencho, glad to hear it worked out. (that google link doesn't work at all).

@Daniel, if you still want the paraphrase: the basic model (Hull example), models the equity as a call option -- with the asset value [$]V_t[$] playing the role of the underlying and a call strike equal to the debt value L. In that model, the company goes bankrupt if and only if the asset value is less than the debt at expiration T. But paths where the intermediate asset values (prior to expiration) fall below the strike (=debt) are fine (do not result in bankruptcy) as long as [$]V_T > L[$]. 

The second model tightens up the bankruptcy criterion so that the company is bankrupt if the asset value i) falls below the barrier H at *any* time prior to T  or ii) is below L at expiration. The natural choice is H = L, but you could argue that H < L could make sense, too. In any event, If H is much smaller than L (H << L), the likelihood of the particle reaching the barrier is small, so it's as if the barrier does not exist, and so you should recover the basic model results.  

Clear. Thanks.
There is a PDE to describe first exit time. For Merton this is OK. For double barriers it should also be feasible? The company's equity would be the derivative.
http://www.datasimfinancial.com

I tell the kids, somebody’s gotta win, somebody’s gotta lose. Just don’t fight about it. Just try to get better.
 
User avatar
Alan
Posts: 8984
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

Re: System of nonlinear equations

June 20th, 2017, 2:30 pm

Well, you'd have to have some plausible boundary condition at the putative upper barrier U. Maybe you could say the company gets acquired there, boosting the equity value to some  value [$]A[$], so [$]E(t,V) = A[$] if [$]V = U[$]. I'd want [$]A \ge U - L[$], where [$]L[$] is the debt; otherwise the problem sounds wonky. (So, this should keep the equity value (= call price) monotone increasing in the asset value [$]V[$] (= underlying price) at a fixed time slice, which might be needed for a solution to the nonlinear system to exist). 
ABOUT WILMOTT

PW by JB

Wilmott.com has been "Serving the Quantitative Finance Community" since 2001. Continued...


JOBS BOARD

JOBS BOARD

Looking for a quant job, risk, algo trading,...? Browse jobs here...