multivariate newton's method not working

Hi. So I'm trying to get something that will tell me the right number of option contracts per dollar owned, for call options of various strike prices, given their current prices per option contract, the strike prices, and the estimated mean and standard deviation of the underlying stock price at expiration. Supposedly Newton's method works in multiple dimensions just as it does in one dimension. That is to say, if you want to close in on a spot where the derivative is 0 of a function which will probably be a minimum or a maximum, you would start from an initial guess x, and then for each additional iteration, Δx = - f'/f'', where f' is the first derivative and f'' is the 2nd derivative of the function you're trying to maximize. In multivariate form, Δx = - (∇²f)^(-1) * ∇f where ∇f is the gradient of the function and ∇²f is the Hessian, that is, the higher dimensional equivalent of the 2nd derivative, which is a matrix of partial derivatives, and (∇²f)^(-1) is the matrix inverse of ∇²f. And here is the code I have created. It doesn't produce any errors, but I don't understand why it's not converging. It seems to be slowly going in the wrong direction in fact. The current settings are a stock that's currently 445 dollars a share and expected to rise $1.50 per day for the next 9 days (cough cough Eli Lily cough cough) but with a standard deviation of $5.30 per day. And standard deviation increases with the square root of time and mean with the time, so in other words after 9 days, the standard deviation of the price's probability distribution is then 3 times as much or 15.90 and the mean price per share expected would be up 13.50 from the current price of 445.

--------------------

clear all
close all

currentprice=445;
optionstrikeprices=[442.5, 460, 470, 480];
calloptionprices=[8.675, 1.88, .685, .4];
expectedpriceatexpiration=currentprice+1.5*9;
standarddeviationatexpiration=5.3*sqrt(9);

%The cumulative probability function of a gaussian random variable
%with standard deviation 1 is .5*(1+erf(x/sqrt(2)))

quantitytoown=zeros(size(optionstrikeprices));
currentutility=0;
utilitydifferentials=zeros(1,size(optionstrikeprices,2)-1);
utilityhessian=zeros(size(utilitydifferentials,2),size(utilitydifferentials,2));
gaussiandistributiondomain=-4:.1:4;
gaussiandistribution=exp(-.5*gaussiandistributiondomain.^2);
gaussiandistribution=gaussiandistribution/sum(gaussiandistribution);

for counter3=1:9
    for counter2=1:size(utilitydifferentials,2)
        quantitytoowntest=quantitytoown+[zeros(1,counter2-1), .00001, zeros(1,size(quantitytoown,2)-counter2-1), -.00001];
        utility=0;
        for counter=1:size(gaussiandistribution,2)
            endingprice=expectedpriceatexpiration+gaussiandistributiondomain(counter)*standarddeviationatexpiration;
            calloptionvaluesatend=endingprice-optionstrikeprices;
            calloptionvaluesatend=calloptionvaluesatend.*(calloptionvaluesatend>0);
            performance=(1+quantitytoowntest*(calloptionvaluesatend-calloptionprices)')*2/(1+endingprice/currentprice); %*2/(1+endingprice/currentprice) because it is making a benchmark out of 50% dollars and 50% shares at their starting price
            if performance<1e-20
                performance=1e-20;
            end
            utility=utility+gaussiandistribution(counter)*log(performance);
        end
        utilitydifferentials(counter2)=(utility-currentutility)*100000; %Now stores the gradient
        for counter4=1:counter2 %Now need the hessian
            quantitytoowntest2=quantitytoowntest+[zeros(1,counter4-1), .00001, zeros(1,size(quantitytoown,2)-counter4-1), -.00001];
            utility2=0;
            for counter=1:size(gaussiandistribution,2)
                endingprice=expectedpriceatexpiration+gaussiandistributiondomain(counter)*standarddeviationatexpiration;
                calloptionvaluesatend=endingprice-optionstrikeprices;
                calloptionvaluesatend=calloptionvaluesatend.*(calloptionvaluesatend>0);
                performance=(1+quantitytoowntest2*(calloptionvaluesatend-calloptionprices)')*2/(1+endingprice/currentprice);
                if performance<1e-20
                    performance=1e-20;
                end
                utility2=utility2+gaussiandistribution(counter)*log(performance);
            end
            utilityhessian(counter2,counter4)=((utility2-utility)*100000-utilitydifferentials(counter2))*100000;
        end
    end
    for counter2=1:(size(utilitydifferentials,2)-1)
        utilityhessian(counter2, (counter2+1):size(utilitydifferentials,2))=utilityhessian((counter2+1):size(utilitydifferentials,2), counter2)';
    end
    utilitydifferentials=utilitydifferentials*utilityhessian^-1;
    quantitytoown=quantitytoown-[utilitydifferentials, -sum(utilitydifferentials)]
    currentutility=0;
    for counter=1:size(gaussiandistribution,2)
        endingprice=expectedpriceatexpiration+gaussiandistributiondomain(counter)*standarddeviationatexpiration;
        calloptionvaluesatend=endingprice-optionstrikeprices;
        calloptionvaluesatend=calloptionvaluesatend.*(calloptionvaluesatend>0);
        performance=(1+quantitytoown*(calloptionvaluesatend-calloptionprices)')*2/(1+endingprice/currentprice);
        if performance<1e-20
            performance=1e-20;
        end
        currentutility=currentutility+gaussiandistribution(counter)*log(performance);
    end
    currentutility
end

quantitytoown = 1×4

1.0e-05 * 0.9971 -0.0000 -0.0027 -0.9944

currentutility = -0.0148

quantitytoown = 1×4

1.0e-04 * 0.0822 -0.0000 -0.1178 0.0356

currentutility = -0.0148

quantitytoown = 1×4

1.0e-04 * 0.0646 0.0000 -0.2354 0.1707

currentutility = -0.0149

quantitytoown = 1×4

1.0e-04 * 0.0471 0.0000 -0.3529 0.3058

currentutility = -0.0149

quantitytoown = 1×4

1.0e-04 * 0.0296 0.0000 -0.4705 0.4409

currentutility = -0.0149

quantitytoown = 1×4

1.0e-04 * 0.0120 0.0000 -0.5880 0.5760

currentutility = -0.0150

quantitytoown = 1×4

1.0e-04 * -0.0055 0.0000 -0.7056 0.7111

currentutility = -0.0150

quantitytoown = 1×4

1.0e-04 * -0.0231 0.0000 -0.8231 0.8462

currentutility = -0.0150

quantitytoown = 1×4

1.0e-04 * -0.0406 0.0000 -0.9407 0.9812

currentutility = -0.0151

[John D'Errico]

Why not use a provided optimization tool, instead of writing your own? It is NEVER a good idea to write your own basic numerical methods code. Unless of course, you understand them so well that you can write it bettter than professionals in the field of numerical analysis?

[Torsten]

I don't see what x and f from Newton's formula are in your code and where x is updated. So it's impossible to gives advice.

[Conrad]

Re: Torsten

Granted I should have maybe commented things more. Ok, so f is "utility". What is utility? It represents the integral, weighted by a gaussian distribution (defined as a vector before the loop starts, named gaussiandistribution), of the log of the effect of the gamble on my total bankroll. This is slightly more useful than just using "expected value" because this will allow one to resolve the question of how MUCH one should bet on a given gamble. If you want to maximize your statistical expectation, you bet as much as you can. The problem is then you lose all your money as soon as things don't go your way. So if you take the LOG of the amount of money you have at the end before finding the expected value, you get something more useful, because the log of 0 is negative infinity, so a complete loss represents something that absolutely must be avoided in the optimization process. Because really one's progress in gambling is measured on a log scale, and if you're betting too big, even if your steps forward are more frequent, if your bets are too big for the money you have, the steps back will be bigger on a log scale.

x on the other hand, represents the number of positions owned of each of the types of options. There are 4 kinds of options but there are only 3 free variables defining it, however, because I'm constraining there to be a short option for every long option, so that the sum of the number of options is 0, and so the first variable is actually long one of the first of 4 options, and short the last one of the 4, the second variable is long the 2nd of the 4 options and short the last one, and the 3rd variable is long the 3rd of the 4 options and short the last one. This is why in calculating the differential, you see something like

quantitytoowntest2=quantitytoowntest+[zeros(1,counter4-1), .00001, zeros(1,size(quantitytoown,2)-counter4-1), -.00001];

This is finding the differential change from adding .00001 to one of the first 4 variables of the vector defining the number of each kind of option to own, and subtracting .00001 to the last one, and then the utility at this new value is found, and the old value at the previous position is subtracted from it and multiplied by 100000, to approximate taking a derivative.

Re: John D'Errico

You seem to be under the impression that it is far more extraordinary to thoroughly understand Newton's method, which I learned in high school calculus class, than it is to understand the use of obscure matlab toolboxes, which I have not learned anywhere and have no way of knowing. The function I'm dealing with is the outcome of a computation with several steps and I don't know anything about this thing you're talking about, other than that someone else wrote it and I can't even imagine how it could "optimize" my thing, because how could IT run my operation to produce the result inside of this other program someone else wrote, it's not a simple function with a simple formula.

What's more, I did something VERY similar with options about 10 years ago so I know multivariate Newton's method SHOULD work for this, because that thing converged very thoroughly in about 10 iterations. So it should be working. Unfortunately that hard drive died with no warning (curse you seagate).

Answer 1

[Matt J]

You seem to be under the impression that it is far more extraordinary to thoroughly understand Newton's method, which I learned in high school calculus class, than it is to understand the use of obscure matlab toolboxes

It is not that it is extraordinary to understand plain-vanilla Newton's method. It is just that plain Newton's method is rarely ever used, as it is less reliable than certain refined versions of the algorithm, some of which are implemented in the Optimization Toolbox for your convenience. In particular, damped versions of Newton's method, which adjust the stepsize adaptively using a line search, will have a larger radius of converge than plain-vanilla Newton.

With that in mind, your problem could be that you are simply not starting close enough to the desired solution for convergence to occur. Like the others have said though, the code is very hard to read, and so any bugs present in your implementation are very likely to escape our notice.

[Conrad]

Well thanks for taking a look, especially Torsten. Damn but LLY really did go up an average of 1.5 dollars per share over those 9 days, almost exactly, right in the center of the bellcurve of my expectations. If only I had figured out the issues with the code and bought its recommendation, I would have really cleaned up. Also I've been hurting all this time with a severe neck injury so I haven't been able to do much in this time anyway, but it looks like no one knows what the problem is any more than I do. I can say however, that I don't think it actually needs to be anywhere near the solution for it to converge, based on the program I wrote 10 years that did much the same thing. I believe the function doesn't have any tricky saddle points in it so it should converge from pretty much anywhere, which is why the initial guess was owning no position. It should have converged pretty quickly, I must have just done something wrong with it that I'm not seeing.

[Matt J]

I believe the function doesn't have any tricky saddle points in it so it should converge from pretty much anywhere

That doesn't matter. Consider the strictly convex function f(x)=|x|^1.5, which has no saddle points. If you initialize plain vanilla Newton's method from , it will just oscillate forever between +1 and -1, never converging to the minimum at x=0.

syms x real

f=abs(x)^(3/2);

nextPoint(x) = simplify( x - diff(f,1)./diff(f,2) ,3)  %newton update step

nextPoint(x) =