SVD routine error?

[wdgann2002]

Hello,

 

I'm new to quantLib but I'm getting some bad numbers back from the SVD getV() routine.  The other matricies are ok but not not V.    Are there any know problems with that code or is the problem (more likely) somewhere on my side?  I did a search through the archive but could not find any  Here is the code segment I'm looking at.  Thanks ahead of time for anyone's help.

 

Rick

 

  Matrix U(A.rows(), A.rows());
  Matrix V(A.columns(), A.columns());
  Matrix S(A.rows(),A.rows());
  Array s(A.columns());
  SVD svd(A);
  svd.getU(U);
  svd.getV(V);      <---------- bad values (they are there but they are wong????)
  svd.getS(S);
  Matrix Vt=transpose(V);
  Matrix Lt = U*S;

---

Do you Yahoo!?
Yahoo! Finance: Get your refund fast by filing online

[Luigi Ballabio-2]

On 2004.02.18 02:16, wdgann2002 wrote:
> I'm new to quantLib but I'm getting some bad numbers back from the  
> SVD getV() routine.  The other matricies are ok but not not V.    Are
> there any know problems with that code or is the problem (more  
> likely) somewhere on my side?

Rick,
        it might very well be that the problem is in the code--it was  
contributed only recently and did not have a lot of testing.

Can you file this report in the bug tracker? This way it doesn't get  
lost in some mailbox in case that we're not able to tackle it right  
away (which is likely.)

Thanks,
        Luigi

[Neil P Firth]

Hi,

The SVD class was my doing. One of my early attempts ;) If you have any
questions about it feel free to email me. If you have some test cases or
other code it would be good to get them into a test class in the
test-suite (on my todo list..), send them to me and I'll see what I can
do.

Cheers,
Neil

---------------------------------------------------
  Neil Firth
  Brasenose College Oxford OX1 4AJ United Kingdom
  Office: 01865 280616
  [hidden email]
  http://www.maths.ox.ac.uk/~firth
---------------------------------------------------
On Wed, 18 Feb 2004, Luigi Ballabio wrote:

> On 2004.02.18 02:16, wdgann2002 wrote:
> > I'm new to quantLib but I'm getting some bad numbers back from the
> > SVD getV() routine.  The other matricies are ok but not not V.    Are
> > there any know problems with that code or is the problem (more
> > likely) somewhere on my side?
>
> Rick,
> it might very well be that the problem is in the code--it was
> contributed only recently and did not have a lot of testing.
>
> Can you file this report in the bug tracker? This way it doesn't get
> lost in some mailbox in case that we're not able to tackle it right
> away (which is likely.)
>
> Thanks,
> Luigi
>
>
> -------------------------------------------------------
> SF.Net is sponsored by: Speed Start Your Linux Apps Now.
> Build and deploy apps & Web services for Linux with
> a free DVD software kit from IBM. Click Now!
> http://ads.osdn.com/?ad_id=1356&alloc_id=3438&op=click
> _______________________________________________
> Quantlib-users mailing list
> [hidden email]
> https://lists.sourceforge.net/lists/listinfo/quantlib-users
>

[Luigi Ballabio-2]

On 2004.02.20 13:07, Neil P Firth wrote:
> The SVD class was my doing. One of my early attempts ;)
> If you have some test cases or other code it would be good to get  
> them into a test class in the test-suite (on my todo list..), send  
> them to me and I'll see what I can do.

Neil,
        I exported the SVD class to Python for ease of test. The  
interface is a bit different, but the calculations come from your  
class. Here's a transcript of a Python session:

>>> from QuantLib import *
>>> A = Matrix([[1,2,3,4],
...            [2,0,2,1],
...            [0,1,0,0]])
>>> svd = SVD(A)

U should be orthogonal:
>>> U = svd.getU()
>>> print U*transpose(U)

1.169252,0.021675,-0.000000
0.021675,1.002776,-0.000000
-0.000000,-0.000000,1.000000

Close, but no cake. The error is probably too high for numerical  
applications.
V should be orthogonal too:

>>> V = svd.getV()
>>> print V*transpose(V)

1.000000,-0.000000,0.000000,0.000000
-0.000000,0.469295,-0.228611,0.443615
0.000000,-0.228611,0.901521,0.191096
0.000000,0.443615,0.191096,1.629184

Ouch. No good.
Now for the final rush:

>>> S = svd.getS()
>>> print U*S*transpose(V)

1.000000,2.027216,3.011724,4.014609
2.000000,-0.013608,1.994138,0.992695
-0.000000,0.469295,-0.228611,-0.284880

Again, close but not enough.

When I used transpose(A) to initialize the SVD, U was bogus, V was good  
(V*transpose(V) == I up to 6 decimals) and U*S*transpose(V) was just as  
good.

For a square matrix (A minus the rightmost column, or A plus another  
row) the results were uniformly good (up to the 6 decimals I see when  
printing.)


May I suggest that the test case for the test-suite might bubble  
upwards a bit in your todo-list? :)
The test should:
a) create a SVD for a matrix A;
b) retrieve U and verify that U*transpose(U) == I up to some tolerance  
(1e-6? 1e-8? For each element? For the L2 norm of the matrix?)
c) retrieve V and do the same
d) retrieve S and while you're at it, verify that it is diagonal and  
its elements equal the singular values
e) verify that U*S*transpose(V) == A up to some tolerance
f) repeat for a few matrices (square, rectangular either way, with  
elements of the same magnitude, with some elements much larger or much  
smaller.)

Thoughts?

Later,
        Luigi

[Bill Shortall]

----- Original Message -----
From: "Luigi Ballabio" <[hidden email]>
To: "wdgann2002" <[hidden email]>
Cc: <[hidden email]>
Sent: Wednesday, February 18, 2004 3:33 AM
Subject: Re: [Quantlib-users] SVD routine error?

> On 2004.02.18 02:16, wdgann2002 wrote:
> > I'm new to quantLib but I'm getting some bad numbers back from the
> > SVD getV() routine.  The other matricies are ok but not not V.    Are
> > there any know problems with that code or is the problem (more
> > likely) somewhere on my side?
>
> Rick,
> it might very well be that the problem is in the code--it was
> contributed only recently and did not have a lot of testing.
>
> Can you file this report in the bug tracker? This way it doesn't get
> lost in some mailbox in case that we're not able to tackle it right
> away (which is likely.)
>
> Thanks,
> Luigi

Hi Luigi,

,  There are some more problems with SVD
In the line

T scale = QL_MAX(QL_MAX(QL_MAX(QL_MAX(
               QL_FABS(s[p-1]),QL_FABS(s[p-2])),QL_FABS(e[p-2])),
               QL_FABS(s[k])),QL_FABS(e[k]));

The Microsoft compiler VC6 doesn't like the macro expansion. Won't work
in release mode. I think its the recursion.
Try something like
   double scale = max_funct(max_funct(max_funct(max_funct(
                       abs_func(s[p-1]),abs_func(s[p-2])),abs_func(e[p-2])),
                       abs_func(s[k])),abs_func(e[k]));

where
   double  max_funct(double a, double b){
 if(a>b) return a;
 else
 return b;
}
and also
double abs_func(double x) // returns |x| {
 if( x <= 0 )
  return -x ;
 else
  return x;
}

This class is going to be a real bear.

Regards..  Bill Shortall

>
>
> -------------------------------------------------------
> SF.Net is sponsored by: Speed Start Your Linux Apps Now.
> Build and deploy apps & Web services for Linux with
> a free DVD software kit from IBM. Click Now!
> http://ads.osdn.com/?ad_id=1356&alloc_id=3438&op=click
> _______________________________________________
> Quantlib-users mailing list
> [hidden email]
> https://lists.sourceforge.net/lists/listinfo/quantlib-users
>

[Neil P Firth]

Hi,

Sorry about the trouble this class is causing ;)

I've checked in a test to testsuite\matrices.cpp which looks at a couple
of cases. The accuracy is good to about 1e-6, which isn't that great.
Anyone having difficulty please contribute your Matrices as test cases.

The algorithm is ported from the TNT / JAMA project, and should have m>=n
(rows >= columns). I don't know why this is however. I was considering
throwing an exception if n>m and forcing the user to transpose the matrix
first.

Are there any design thoughts on extending the linear algebra. For
instance, I'd like to be able to swap the TNT SVD class for a NAG or
CLAPACK routine to utilise the BLAS and the shared memory processors on
the Oxford computational finance cluster.

As for the macro problem:

>
> ,  There are some more problems with SVD
> In the line
>
> T scale = QL_MAX(QL_MAX(QL_MAX(QL_MAX(
>                QL_FABS(s[p-1]),QL_FABS(s[p-2])),QL_FABS(e[p-2])),
>                QL_FABS(s[k])),QL_FABS(e[k]));
>
> The Microsoft compiler VC6 doesn't like the macro expansion. Won't work
> in release mode. I think its the recursion.
> Try something like
>    double scale = max_funct(max_funct(max_funct(max_funct(
>                        abs_func(s[p-1]),abs_func(s[p-2])),abs_func(e[p-2])),
>                        abs_func(s[k])),abs_func(e[k]));

I think it would be best to do half the work in
double temp_scale =....
then do the rest in
double scale =...
as we want to use the QL_MAX macros everywhere for maintainability.

Neil

[Luigi Ballabio-2]

Hi Neil,

At 23:39 22/02/2004, Neil P Firth wrote:
>Sorry about the trouble this class is causing ;)

No problem.

>I've checked in a test to testsuite\matrices.cpp which looks at a couple
>of cases. The accuracy is good to about 1e-6, which isn't that great.
>Anyone having difficulty please contribute your Matrices as test cases.
>
>The algorithm is ported from the TNT / JAMA project, and should have m>=n
>(rows >= columns). I don't know why this is however. I was considering
>throwing an exception if n>m and forcing the user to transpose the matrix
>first.

I'm having a quick look at the thing. Question: line 122:

             if (wantu & (k < nct)) {

Are you sure it shouldn't be && ?


>Are there any design thoughts on extending the linear algebra. For
>instance, I'd like to be able to swap the TNT SVD class for a NAG or
>CLAPACK routine to utilise the BLAS and the shared memory processors on
>the Oxford computational finance cluster.

Ehm, I'm not exactly sure what you mean. "Extending", as in ...?

>As for the macro problem:
> >
> > ,  There are some more problems with SVD
> > In the line
> >
> > T scale = QL_MAX(QL_MAX(QL_MAX(QL_MAX(
> >                QL_FABS(s[p-1]),QL_FABS(s[p-2])),QL_FABS(e[p-2])),
> >                QL_FABS(s[k])),QL_FABS(e[k]));
> >
> > The Microsoft compiler VC6 doesn't like the macro expansion. Won't work
> > in release mode. I think its the recursion.

This is weird. I've been using VC6 for a long time and it didn't object to
compiling the above... Bill, what is your configuration?

Later,
         Luigi

[Neil P Firth]

Hi,

>
> I'm having a quick look at the thing. Question: line 122:
>
>              if (wantu & (k < nct)) {
>
> Are you sure it shouldn't be && ?
>
Does look rather like it should. I'll email the guys at the TNT / JAMA
project and see what they say about it. I'll also ask what accuracy to
expect.

>
> >Are there any design thoughts on extending the linear algebra. For
> >instance, I'd like to be able to swap the TNT SVD class for a NAG or
> >CLAPACK routine to utilise the BLAS and the shared memory processors on
> >the Oxford computational finance cluster.
>
> Ehm, I'm not exactly sure what you mean. "Extending", as in ...?
>
I mean that I want to make it easy to swap from using the current SVD
class to a NAG SVD routine, or CLAPACK SVD routine, if they are available.
Would 'pluggable' be a better word..

Neil

[Luigi Ballabio-2]

Hi Neil,

         some news:

>I've checked in a test to testsuite\matrices.cpp which looks at a couple
>of cases. The accuracy is good to about 1e-6, which isn't that great.
>Anyone having difficulty please contribute your Matrices as test cases.

Actually, there was an error in the definition of norm() in the test.
With the correct definition your routine seems to be good to 1e-12.

The bad news is that the error in norm() was hiding a bug in pseudoSqrt().
Apparently, the thing couldn't calculate a sqrt to save its life.
Nando, can you take a look at it?


>The algorithm is ported from the TNT / JAMA project, and should have m>=n
>(rows >= columns). I don't know why this is however. I was considering
>throwing an exception if n>m and forcing the user to transpose the matrix
>first.

No, I just extended the thing. You just decompose the transposed matrix,
and swap the results. The proof is in a comment in svd.cpp. You might want
to have a look and check that it's ok. But I transposed one of the matrices
in the test, and it still works.

Oh, and while I was at it, I quantlibified the code a bit more. I hope
you've no problems with this.


> > Ehm, I'm not exactly sure what you mean. "Extending", as in ...?
> >
>I mean that I want to make it easy to swap from using the current SVD
>class to a NAG SVD routine, or CLAPACK SVD routine, if they are available.
>Would 'pluggable' be a better word..

Encapsulating the NAG/CLAPACK routine in a class with the same interface as
the current one would go a long way towards your goal. Once it is done, one
can just pass the class to use as a template argument.

Later,
         Luigi

[Ferdinando M. Ametrano-3]

Luigi wrote:
>Actually, there was an error in the definition of norm() in the test.
>With the correct definition your routine seems to be good to 1e-12.
>
>The bad news is that the error in norm() was hiding a bug in pseudoSqrt().
>Apparently, the thing couldn't calculate a sqrt to save its life.
>Nando, can you take a look at it?
it's now fixed in the CVS: the error in norm() was hiding a bug in the test
case code, not in pseudoSqrt()

ciao -- Nando

[daniel saitowitz]

Is there any code for double barrier options?
I can't seem to find any, if not maybe I can contribute some.

Regards
Daniel

[Luigi Ballabio-2]

Hi Daniel,

On 2004.02.26 00:45, daniel saitowitz wrote:
> Is there any code for double barrier options?

No, there's not.

> I can't seem to find any, if not maybe I can contribute some.

That would be great.

Thanks,
        Luigi

[daniel saitowitz]

Hi Guys,

I have prototyped some methods for Double Barriers in VBA and now
looking at doing a QuantLib integration.

I propose a design similar to the Current Single Barriers:

Instrument:
        doublebarrieroption.hpp
        doublebarrieroption.cpp

Pricing Engine:
        analyticdoublebarrierengine.hpp
        analyticdoublebarrierengine.cpp
        trinomialdoublebarrierengine.hpp
        trinomialdoublebarrierengine.cpp

Test Suite:
        doublebarrieroption.hpp
        doublebarrieroption.cpp


Implementations are from Haug as a series of Single Barriers so can
reuse some of the single barrier stuff here.
For the Tian's trinomial method, I see there is a trinomialtree class,
and not quite sure how to integrate with it. The tree for double
barriers has different probabibilities for top half , middle and bottom.
So can I implement it straight in the pricing engine or is it necessary
to extend the trinomialtree class?

I don't see an implementation of Tian's trinomial for Single Barriers,
so maybe I can also add:
        trinomialbarrierengine
Plus some optimisations using Richardson extrapolation. Also on the
analytic side will implement the BGK barrier adjustments.

I saw some nice uml on the website, what is being used to model this?
Any other hints or tips on getting started appreciated.

Regards
Daniel

[Ferdinando Ametrano-3]

Hi Daniel

>I have prototyped some methods for Double Barriers in VBA and now
>looking at doing a QuantLib integration.
>
>I propose a design similar to the Current Single Barriers:
>[...]
ok

>Implementations are from Haug as a series of Single Barriers so can
>reuse some of the single barrier stuff here.
please include as many numerical examples from Haug's book and
VBA/spreadsheet as possible in the test suite, plus logical checks as
knock-in+knock-out=knock-less, etc.

>For the Tian's trinomial method, I see there is a trinomialtree class,
>and not quite sure how to integrate with it.
the only trinomial tree available in QuantLib it's the one used for
single-factor interest rate model. If I'm not wrong it's the Hull-White
trinomial tree.
The Tian tree already available in QuantLib is the binomial one (third
moment matching)

>  The tree for double
>barriers has different probabibilities for top half , middle and bottom.
>So can I implement it straight in the pricing engine or is it necessary
>to extend the trinomialtree class?
I'm not familiar with the problem, but I guess the tree could be reused in
other pricing engines or replaced by different trees in the double barrier
engine. If this is the case it would be better not to have its
implementation in the pricing engine: see how binomial trees are used in
the binomial pricing engine

>I don't see an implementation of Tian's trinomial for Single Barriers,
>so maybe I can also add:
>         trinomialbarrierengine
sure!

>Plus some optimisations using Richardson extrapolation.
I look forward to your proposed design as I'm interested to enable
Richardson extrapolation in the (future) Finite Difference engines.


>  Also on the analytic side will implement the BGK barrier adjustments.
you're welcome.

Please pay special attention to the test-suite: I am really worried about
contributions without an appropriate test-suite. Please try to reference
the relevant papers/books and to reproduced all numerical examples provided
in your sources

>I saw some nice uml on the website, what is being used to model this?
I'm sure Luigi will step in (at least) here...

ciao -- Nando

[Luigi Ballabio-2]

Hi,

>> I saw some nice uml on the website, what is being used to model  
>> this?
> I'm sure Luigi will step in (at least) here...

I've been using Dia (http://www.lysator.liu.se/~alla/dia/)
If I were to need them for a paper instead of a web page, I'd probably  
use the pst-uml package for LaTeX (you can Google for it).

Later,
        Luigi

[daniel saitowitz]

>Please pay special attention to the test-suite: I am really worried
about
>contributions without an appropriate test-suite. Please try to
reference
>the relevant papers/books and to reproduced all numerical examples
provided
>in your sources

Sure - I totally understand the importance of this!
Both test-suite and references.
I have good data from the Tian paper, unfortunately the Haug paper and
Luo paper, don't have enough data for replication, so I am benchmarking
them against Tian and some other papers.


Regards
Daniel

[rjcookusa]

Haug type, analytic double barrier options

Hi

I created a double barrier class based on the barrier option class some time ago, but I don't think it ever made it into CVS. The implementation is very much like the HAUG VBA model, but I think I took out curved barriers etc. If you need to look at it, I can find it.

Julian Cook

[daniel saitowitz]

Hi Julian,

I would appreciate that. If you can dig it out and send it - I would be
interested.

Regards
Daniel

-----Original Message-----
From: [hidden email]
[mailto:[hidden email]] On Behalf Of
[hidden email]
Sent: 02 March 2004 22:32
To: [hidden email]
Subject: [Quantlib-users] RE: Double Barrier Options


Haug type, analytic double barrier options

Hi

I created a double barrier class based on the barrier option class some
time ago, but I don't think it ever made it into CVS. The implementation
is very much like the HAUG VBA model, but I think I took out curved
barriers etc. If you need to look at it, I can find it.

Julian Cook


-------------------------------------------------------
SF.Net is sponsored by: Speed Start Your Linux Apps Now.
Build and deploy apps & Web services for Linux with
a free DVD software kit from IBM. Click Now!
http://ads.osdn.com/?ad_id=1356&alloc_id=3438&op=click
_______________________________________________
Quantlib-users mailing list [hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-users

[rjcookusa]

This is the info I sent to Daniel Saitowitz without the source code
attachments for QuantLib.Pricers.DblBarrierOption::DblBarrierOption

Daniel:
These are the source files (dblbarrieroption.cpp etc). They were last
compiled into the QuantLib-0.3.0 source tree. There is also a reference to the model
in quantlib.hpp, so I attached that.

I also created a test case that goes in the examples directory i.e. I created
a new dir called D:\QuantLib-0.3.0\Examples\BarrierOption and created a VC++
dsw for it. I ve attached the actual test file (BarrierOption.cpp) that was
modified from another test case. If you want the entire directory I can send it
as well. The test case actually diffs the value with the VBA prototype, which
I previously verified as correct vs Fenics.

It should look like this

D:\QuantLib-0.3.0\Examples\BarrierOption\build\Debug>BarrierOption
Time to maturity = 0.25
Underlying price = 100
Strike = 100
Risk-free interest rate = 0.1
Volatility = 0.15
BarType = 1
LoBarrier = 90
HiBarrier = 110

Method              Value                       EstimatedError
Black Scholes   2.91752076761027    0.00000000000000    0.0000  0.0000
Double Barrier  1.20547906151793    0.00000010708526    0.0000  0.0000

It should be noted that the VBA price = 1.20547916860319, which is not
exactly the same as the C++ price here.

regards

Julian Cook