4

I have problems when I'm trying to do numerical differentiation in Matlab. But my question might me more about numerical analysis than about Matlab.

I have an array with 9 data points that represents f(x) for 9 different x. I need to find f''(x) numerically. The values I have for x and f(x) are

x = [2271.38, 2555.30, 2697.26, 2768.24, 2839.22, 2910.20, 2981.18, 3123.14, 3407.06]

f(x) = [577.4063, 311.3341, 193.0833, 141.3048, 95.1501, 58.8130 32.4931, 6.9511, 0.1481]

and I can interpolate to get a smooth curve. I use spline interpolation but is some other interpolation preferable when you are going to differentiate?

I've tried different methods:

Just simple forward,backward and central difference quotients

A wavelet based method: http://www.mathworks.com/matlabcentral/fileexchange/13948-numerical-differentiation-based-on-wavelet-transforms

and the derivest suite: http://www.mathworks.com/matlabcentral/fileexchange/13490-adaptive-robust-numerical-differentiation

Non of this have worked satisfactory. The second derivative is very unstable with respect to the step length and the adaptive method in the derivest suite works terribly bad. Maybe I'm just using it in the wrong way!

Any help is appreciated!

Thanks in advance

DoubleTrouble
  • 902
  • 2
  • 7
  • 19
  • 1
    What are those data points ? Are they actual measurements with errors, are they test data to check data points ? The function itself looks like exponential decay. – Thorsten S. Mar 28 '12 at 11:01
  • Thank you for your reply Thorsten! They are actual measurements with error and here is the interpolated "function": http://i.imgur.com/j3AZl.png – DoubleTrouble Mar 28 '12 at 11:11
  • Given your reply under High Performance Marks answer "They are just data points from a stock market" my crystal ball says that your supervisor is trying to forecast stock movements and you are trying to get the derivative at the last point which fails with each method you have been used so far...? – Thorsten S. Mar 28 '12 at 12:27
  • 2
    No, good guess though! I'm actually trying to recover the risk-neutral density from option prices by using a well known relation by Dupire namely that the density function is equal to the discounted second derivative of call option prices! – DoubleTrouble Mar 28 '12 at 14:49

2 Answers2

9

I imagine it was you who asked a similar question on MATLAB Central the other day. You did not post your data there, so I could not really give a good answer then.

Estimation of the second derivative is a difficult thing to do. It is an ill-posed problem. Differentiation is itself a noise amplifier, so estimating the second derivative is "twice" as bad. It is simply not an easy thing to do, certainly not well.

Using this set of points, I chose to estimate a spline model using my SLM toolbox.

x = [2271.38, 2555.30, 2697.26, 2768.24, 2839.22, 2910.20, 2981.18, 3123.14, 3407.06];
f = [577.4063, 311.3341, 193.0833, 141.3048, 95.1501, 58.8130 32.4931, 6.9511, 0.1481];

First of all, plot the data. What can I learn from that plot? Are there any inferences I might choose to make?

pure data plot

The simple plot tells me, along with your comments, that I expect this function to be a monotonic decreasing function. It appears to be asymptotically linear at each end, like a hyperbolic segment, with positive curvature over the entire domain.

So now I would use this information to build a model for your data, using my SLM toolbox.

slm = slmengine(x,f,'plot','on','decreasing','on','knots',20, ...
    'concaveup','on','endconditions','natural');

slmengine is designed to take information from you, in the form of prescriptions for the shape of the curve. What you will find is that by providing such information, it greatly regularizes the shape of the result, to conform to your knowledge of the process. Here I was just making a few guesses about the curve shape from your comments.

In the above call, I instructed SLM to:

  • make a plot of the result when done
  • create a monotone decreasing function of x
  • use 20 equally spaced knots
  • force the curve to have an everywhere positive second derivative
  • set the second derivatives at the end to be zero

The plot as generated is a gui itself, allowing you to plot the function and data, but also to plot derivatives of the result. The vertical green lines are the knot locations.

curve and data

Here we see that the curve fit is a reasonable approximation to what you are looking for.

How about the second derivative plot? Of course, SLM is a piecewise cubic tool. Therefore the second derivatives are only piecewise linear. Is this a problem? Will you ask me to provide a tool for higher order splines? Too bad, but no, I won't. Those higher order derivatives are far too poorly estimated to ask for a highly smooth result. In fact, I would be quite happy with this prediction. Note that the glitches in the second derivative were consistent. If I used more knots or fewer, they were still there. This is a nice way to know if the shape you see is a feature of the curve, or merely an artifact of knot placement.

See that the constraints I placed on the shape of the curve resulted in a fit that was quite reasonable, despite the fact that I used far more knots than I had data points. SLM had no trouble in the estimation.

2nd derivative

If I wish to try for a smoother estimate of the second derivative, simply use more knots. SLM is relatively fast. So with 50 knots, we get a very similar result for the second derivative curve.

2nd derivative, 50 knots

You can find SLM (here) on MATLAB Central. It does require the optimization toolbox.

  • Oh my god John! I have now spent an hour with your SLM tool and it truly mops the floor with all other tools! I got the result I expected and I thank you dearly for taking your time to write this extensive reply to my question. P.S. It was probably not me who asked the other day. I did however ask this question at Matlab Central today! – DoubleTrouble Mar 28 '12 at 14:45
  • Very nice answer with great explanation of not only how to do this, but also how to interpret the result. +1 since I can't give +5. – Jonas Mar 28 '12 at 16:55
1

This is an extended comment, not an answer:

What do you know about the function that generates these points ? If, for example, you have good reason to believe that it is a degree-2 polynomial then your first step would be to find the degree-2 polynomial that fits your points best, then take the 2nd derivative of it. If you think that a spline (of some sort) is a better description of f, then fit a spline.

I haven't, and I'm not going to, plotted the graph of your function so I'm not going to hazard a guess as to what would be a reasonable sort of function to start with.

If you have no knowledge of the function then you can fit whatever curve gives you the results that you like best -- just don't kid yourself that you have got the 2nd derivative of the function. You'll have the 2nd derivative of the function you have chosen but will have no more knowledge about the 2nd derivative of the function that created your data than you had at the start.

High Performance Mark
  • 77,191
  • 7
  • 105
  • 161
  • Thank you for your comment and for making me realize that I left out many details in my original post. Unfortunately, I know nothing about the "function" f. They are just data points from a stock market. I plotted the graph of f(x), here it is: http://i.imgur.com/j3AZl.png The second derivative SHOULD look something like a Gauss bell curve. I get something like that with the derivest suite, but it's too rough! Could that be because I only have 9 data points to start off with? – DoubleTrouble Mar 28 '12 at 11:06