Help: Speed up for loop with two functions

Hello everyone, i am currently running this code. But unfortunately it takes for ever to finsish. A 12x10 matrix took nearly 30 minutes and i need one for 1415596x10. Is there any way i can speed up the process? Is there something that i am missing?

My goal is to optimize my estimation parameters of fun with a particle swarm algorithm. For each day in my sample (6009 days), i have multiple data. I want to minimize the error between observed data (implied_volatility) and my estimation data (via fun) with my 10 parameters.

A = [log_moneyness maturity]; xdata = A; ydata = implied_volatility; 
fun = @(x,xdata) (x(1)+x(2).*xdata(:,2)) + (x(3)+x(4).*xdata(:,2)).*((x(5)+x(6).*xdata(:,2)).*...
    (xdata(:,1) -(x(7)+x(8).*xdata(:,2))) + sqrt((xdata(:,1) - (x(7)+x(8).*xdata(:,2)).^2 ...
    + (x(9)+x(10).*xdata(:,2)).^2))); %original function
parameters_forall = table2array(T_param (:,2:11)); %1415596x10 matrix
%implied_volatility is a 1415596x1 matrix with original values from market
%now i want to minimize the error of my estimate by:   
tDays = size(parameters_forall(:,1));
%%
for i = 1:tDays
    x = parameters_forall(i,:);
    residue = @(x) (sum(fun(x,xdata)-implied_volatility(i)).^2);
    y = particleswarm(residue, 10);
    swarm_parameters(i,:) = y;  
end

Grateful for any help, all the best to you!

[Walter Roberson]

    x = parameters_forall(i,:);
    residue = @(x) (sum(fun(x,xdata)-implied_volatility(i)).^2);

You do not use the x that you assign to there. The x in the next line is a different x.

Also, the way your code is structured, all 1415596x10 of xdata would be used in every optimization, with the only difference between optimizations being the implied volatility. It would seem more likely to me that you would only want xdata to be a subset of the available data, such as the data for one day (possibly less)

[Kai Koslowsky]

So would it make more sense to do it like this:

fun = @(x,xdata) (x(1)+x(2).*xdata(:,2)) + (x(3)+x(4).*xdata(:,2)).*((x(5)+x(6).*xdata(:,2)).*...
        (xdata(:,1) -(x(7)+x(8).*xdata(:,2))) + sqrt((xdata(:,1) - (x(7)+x(8).*xdata(:,2)).^2 ...
        + (x(9)+x(10).*xdata(:,2)).^2)));
%nDays=Number of Days in my sample
%N = Number of data for each day (with N(1)=0)
%xdata as my two input factors to fun
%ydata as my single observed output to fun
%x as all my 10 parameters
for i=1:nDays
    kk = N(i);
    ii = N (i+1);
    Xdata = xdata(kk:ii,:);
    Ydata = ydata(kk:ii,:);
    residue = @(x) (sum(fun(x,Xdata)-Ydata).^2);
    y = particleswarm(residue, 10,[],[],[]);
    swarm_parameters(i,:) = y; 
end

So that the particleswarm minimizes my x parameters, via minimizing the residue/difference in Ydata and fun?

[Walter Roberson]

Yes... except you seem to have lost the implied volatility, which was a scalar ??

[Kai Koslowsky]

%ydata=implied_volatility
for i=1:nDays
    kk = N(i);
    ii = N (i+1);
    Xdata = xdata(kk:ii,:);
    Ydata = ydata(kk:ii,:);
    residue = @(x) sqrt((1/ii).*(sum(fun(x,Xdata)-Ydata).^2));
    y = particleswarm(residue, 10,lb,ub,[]);
    swarm_parameters(i,:) = y; 
end

This is what i worked with now. I had to finalize my residue function. ydata is the implied volatility from obeserved market data, which i calculated before. The parameters, which i get from fun, give me an implied volatility without arbitrage (SVI-method). With particle swarm i try to minimize the estimation error of my parameters. Does this answer your question? :) Or am is still doing something wrong?

[Walter Roberson]

I guess that should be okay.

As discussed before, though, I think your upper bound should be N(i+1)-1 . Otherwise each N(K) except the first and last are used twice, once at the end row for ii = K-1, and once as the start row when ii = K

About how many entries are expected per day?

When I was testing symbolically last night, I found the unexpected conclusion that if you

syms XD [1, size(Xdata,1)]
syms YD [size(Xdata,1), 1]
syms X [1 10]
R = sum(fun(X,XD)-YD).^2))
part1 = solve(diff(R,X(1)),X(1))
R2 = subs(R, X(1), part1)

then the result came out as R2 = 0. Same if you solve the derivative with respect to X(2) .

So at the optimal X(1) or optimal X(2) [according to calculus], the entire residue is 0.

This did not happen for the other variables that I happened to test.

I think the implication is that there might be a lot of different solutions -- which might make it difficult to minimize.

[Kai Koslowsky]

You are absolutely right, about having to adjust N. What i therefore did was:

E = grp2idx(SPX.date); %1415596x1 dates
t = unique (E); %6009 unique dates
Nt = histc(E(:),t); % gives me 6009x1 with number of occourences each day.
null = 0;
N = [null;N];% N now has the size 6010x1
nDays = size(Nt);%6009
for i=1:nDays
    kk = N(i)+1;
    ii = N (i+1);
    Xdata = xdata(kk:ii,:);
    Ydata = ydata(kk:ii,:);
    residue = @(x) sqrt((1/ii).*(sum(fun(x,Xdata)-Ydata).^2));
    y = particleswarm(residue, 10,lb,ub,[]);
    swarm_parameters(i,:) = y; 
end

So instead of subtracting from "ii", i add +1 to "kk". This does the trick as well.

My option data set goes from 1996-2020. In the early years, the average number of options ranges from 10-40 a day, in the later years it increases to 900-2000 per day (each with different maturity and strike price).

I do not however understand the meaning of your second part. When the residue is 0, does that mean it is a optimal estimation parameter?

Maybe a quick background explanation to what i am doing overall. I observed Call/Put data on the SP500 from the market, evaluated implied volatilities of these quotes, Via my function i try to fit an interpolant to them (via SVI Method by Gatheral), with adjustments for maturity and log-moneyness. Now i need to evaluate the interpolant at maturity of 30 calender days and a fine grid of strike prices (No idea yet on how to do so). Then i will interpolated my estimated implied volatilites back to option prices and will try to calculate the risk-neutral density via the work of Breeden/Litzenberger. In the end i will compare the risk neutral PDF to the PDF i observed at the market.

[Walter Roberson]

residue = @(x) sqrt((1/ii).*(sum(fun(x,Xdata)-Ydata).^2));

At the moment, I do not undertand why the residue should be proportional to the index of the last entry? The implication is that if you were to have exactly the same data in (say) January 2018 and (say) March 2018, that the residue should be smaller for the second case because it would have a greater ii and so 1/ii would be smaller ??

On the other hand, the 1/ii factor acts as a constant multiple, and for any two proposed x vectors for the same ii, the relative residue is what is important for determining which is smaller residue, and the relative order is not affected by 1/ii . Likewise, unless you had imaginary components, sum of squared would always be positive, and the relative order of sqrt() of two positive numbers is the same as the relative order of the two numbers itself.

What you need to know for minimization is the relative successes. So you might as well use

residue = @(x) sum((fun(x,Xdata)-Ydata).^2);

Please note that this has an important change to your formula. Watch the location of the () !!

      fun(x,Xdata)           %invoke function to obtain prediction
                - Ydata      %subtract actual value
    (                  ).^2  %square the difference
sum(                       ) %sum the squares

whereas your formula had

        fun(x,XData)             %invoke function to obtain prediciton
                    - Ydata      %subtract actual value
    sum(                    )    %sum the differences
                             .^2 %square the sum

Note: You never store the residues, so you do not care about the absolute values, only about the relative values.

[Kai Koslowsky]

I have attached the formula, which i am trying to replicate. Maybe this will make clear, of what i am triying to accomplish. I agree, that the 1/ii term is misused and i have to adjust it. I will try to understand the logic of what you have explained tomorrow and will get back to this. Many thanks again.

Here Nt is the number of options for each day. The first sigma term are my observed implied volatilites, the second sigma is the estimated implied volatility, which i observe via this function:

where a=a0+a1*tau(time to maturity), same for b,rho,m and sigma, and x=log moneyness of the observed functions. I think, i might have to adjust

%1/ii into 1/(ii-kk)

In general: Would there be any advantage for me, to store each residue for each observation day?

EDIT:

And i think i might have just realized, that i will need to add another sqrt() around function term, as my function calculates the implied variance and not the volatility.

[Walter Roberson]

you currently only store the coefficients, and for coefficients you do not need the sqrt or the division.

You might want to know what the estimate is with the final coefficients. You can evaluate the residue with the coefficients, divide, sqrt, save the value.

ii-kk+1 is the proper divisor

[Kai Koslowsky]

So do i understand you correctly, that the resulting parameters do not change with or without divisor and sqrt? I thought, that when particleswarm is trying to minimize the residue, it will minimize the result via the best parameters and not only the value of the function, so that it fits to ydata. How can i make sure, that my given formula is regarded in each point (sqrt and divisor)?

Edit

For the divisor i agree. I am rather lost on how to make sure, that my parameters are minimized by the function i sent (see Comment above).

[Kai Koslowsky]

Dear Mr. Robinson,

would you like to post something as a answer below, so i can accept your answer? This is the least that i can do. I am truely grateful for your help. I have looked over all your comments many times, and it helped me a lot.

All the best to you!