I have written that this function is part of a research project that involves analyzing time series data from stochastic processes. We have a small number (from 1 to 3) of independent observations of a scalar time series. The observations have different lengths, and each one contains approximately $ 10 ^ 4-10 ^ 5 $ data points The function below `nKBR_moments.m`

it takes an array of observations cells as input, along with other configurations, and generates statistical quantities known as "moments of conditional increments". These are the variables. `M1`

Y `M2`

. For more details of the theory, this research paper describes a similar method.

For research purposes, the function will eventually be evaluated tens of thousands of times, on a desktop computer. An evaluation of this function takes approximately 3 seconds with the test script that I have provided below. Thoughts on optimizing code performance, memory usage or scalability are appreciated.

MATLAB function:

```
function (Xcentre,M1,M2) = nKBR_moments(X,tau_in,Npoints,xLims,h)
%Kernel based moments, n-data
%
% Notes:
% Calculates kernel based moments for a given stochastic time-series.
% Uses Epanechnikov kernel with built in computational advantages. Uses
% Nadaraya-Watson estimator. Calculates moments from n sources of data.
%
%
% Inputs:
% - "X" Observed variables, cell array of data
% - "tau_in" Time-shift indexes
% - "Npoints" Number of evaluation points
% - "xLims" Limits in upper and lower evaluation points
% - "h" Bandwidth
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Processing
dX = (xLims(2)-xLims(1))/(Npoints-1); % Bins increment
Xcentre = xLims(1):dX:xLims(2); % Grid
heff = h*sqrt(5); % Effective bandwidth, for setting up bins
eta = floor(heff/dX+0.5); % Bandwidth for bins optimizing
% Epanechnikov kernel
K= @(u) 0*(u.^2>1)+3/4*(1-u.^2).*(u.^2<=1);
Ks = @(u) K(u/sqrt(5))/sqrt(5); % Silverman's definition of the kernel (Silverman, 1986)
Kh = @(u) Ks(u/h)/h; % Changing bandwidth
% Sort all data into bins
Bextend = dX*(eta+0.5); % Extend bins
edges = xLims(1)-Bextend:dX:xLims(2)+Bextend; % Edges
ndata = numel(X); % Number of data-sets
Xloc = cell(1,ndata); % Preallocate histogram location data
nXdata = cellfun(@numel,X); % Number of x data
key = 1:max(nXdata); % Data key
for nd = 1:ndata
(~,~,Xloc{nd}) = histcounts(X{nd},edges); % Sort
end
Xbinloc = eta+(1:Npoints); % Bin locations
BinBeg = Xbinloc-eta; % Bin beginnings
BinEnd = Xbinloc+eta; % Bin beginnings
% Preallocate
Ntau = numel(tau_in); % Number of time-steps
(M1,M2) = deal(zeros(Ntau,Npoints)); % Moments
(iX,iXkey,XU,Khj,yinc,Khjt) = deal(cell(1,ndata)); % Preallocate increment data
% Pre calculate increments
inc = cell(Ntau,ndata);
for nd = 1:ndata
poss_tau_ind = 1:nXdata(nd); % Possible time-shifts
for tt = 1:Ntau
tau_c = tau_in(tt); % Chosen shift
tau_ind = poss_tau_ind(1+tau_c:end); % Chosen indices
inc{tt,nd} = X{nd}(tau_ind) - X{nd}(tau_ind - tau_c);
end
end
% Loop over evaluation points
for ii = 1:Npoints
% Start and end bins
kBinBeg = BinBeg(ii);
kBinEnd = BinEnd(ii);
% Data and weights
for nd = 1:ndata
iX{nd} = and(kBinBeg<=Xloc{nd},Xloc{nd}<=kBinEnd); % Data in bins
iXkey{nd} = key(iX{nd}); % Data key
XU{nd} = X{nd}(iX{nd}); % Unshifted data
Khj{nd} = Kh(Xcentre(ii)-XU{nd}); % Weights
end
% For each shift
for tt = 1:Ntau
tau_c = tau_in(tt); % Chosen shift
% Get data
for nd = 1:ndata
XUin = iXkey{nd}; % Unshifted data indices
XUin(XUin>nXdata(nd)-tau_c) = (); % Clip overflow
yinc{nd} = inc{tt,nd}(XUin); % Increments
Khjt{nd} = Khj{nd}(1:numel(yinc{nd})); % Clipped weight vector
end
% Concatenate data
ytt = (yinc{:});
Khjtt = (Khjt{:});
% Increments and moments
sumKhjtt = sum(Khjtt);
M1(tt,ii) = sum(Khjtt.*ytt)/sumKhjtt;
y2 = (ytt - M1(tt,ii)).^2; % Squared (with correction)
M2(tt,ii) = sum(Khjtt.*y2)/sumKhjtt;
end
end
end
```

MATLAB test script (no comments are required for this):

```
%% nKBR_testing
clearvars,close all
%% Parameters
% Simulation settings
n_sims = 10; % Number of simulations
dt = 0.001; % Time-step
tend1 = 40; % Time-end, process 1
tend2 = 36; % Time-end, process 1
x0 = 0; % Start position
eta = 0; % Mean
D = 1; % Noise amplitude
gamma = 1; % Drift slope
% Analysis settings
tau_in = 1:60; % Time-shift indexes
Npoints = 50; % Number of evaluation points
xLims = (-1,1); % Limits of evaluation
h = 0.5; % Kernel bandwidth
%% Simulating
t1 = 0:dt:tend1;
t2 = 0:dt:tend2;
% Realize an Ornstein Uhlenbeck process
rng('default')
ex1 = exp(-gamma*t1);
ex2 = exp(-gamma*t2);
x1 = x0*ex1 + eta*(1-ex1) + sqrt(D)*ex1.*cumsum(exp(gamma*t1).*(0,sqrt(2*dt)*randn(1,numel(t1)-1)));
x2 = x0*ex2 + eta*(1-ex2) + sqrt(D)*ex2.*cumsum(exp(gamma*t2).*(0,sqrt(2*dt)*randn(1,numel(t2)-1)));
%% Calculating and timing moments
tic
for ns = 1:n_sims
(~,M1,M2) = nKBR_moments({x1,x2},tau_in,Npoints,xLims,h);
end
nKBR_moments_time = toc;
nKBR_average_time = nKBR_moments_time/n_sims
%% Plotting
figure
hold on,box on
plot(t1,x1)
plot(t2,x2)
xlabel('Time')
ylabel('Amplitude')
title('Two Ornstein-Uhlenbeck processes')
figure
subplot(1,2,1)
box on
plot(dt*tau_in,M1,'k')
xlabel('Time-shift, tau')
title('M^{(1)}')
subplot(1,2,2)
box on
plot(dt*tau_in,M2,'k')
xlabel('Time-shift, tau')
title('M^{(2)}')
```

The test script will create two figures similar to the following.