Lyapunov exponent from time series data for DDE?

If I have 2000 data points for the following equation ,How to find out the Lyapunov exponent for different values of [Tau] .This is a delay differential equation i.e infinite dimensional system.

b = 1;
V0 = 300;
I0 = 0.001;
N1 = 7/4;
p1 = (\[Pi]*b*N1)/V0;
a1 = 1.001; a2 = 0.123; a3 = -3.622*10^-3; b1 = 0.001959; b2 = 0.031; 
b3 = 0.003241; G = 0.5*10^-5; \[CapitalOmega] = 0; C1 = p1/2000;
f = 1*10^3;
Is = 0.0;
\[Tau] = 15;
NL = 1000;

sol1[t_] = 
 NDSolve[{x'[t] - (I0*p1)/(
      2*C1)*((a1*x[t - \[Tau]] - G/(p1*I0)*x[t]) - 
        3/4*x[t - \[Tau]]^3*a2 - 5/8*a3*x[t - \[Tau]]^5) - 
     Is*p1*Cos[y[t]] == 0, 
   y'[t] - \[CapitalOmega] - (I0*p1)/(
      2*C1*x[t])*(b1*x[t - \[Tau]] + 3/4*x[t - \[Tau]]^3*b2 - 
        5/8*b3*x[t - \[Tau]]^5) + (Is*p1)/(2*C1*x[t])*Sin[y[t]] == 0, 
   x[0] == 0.003, y[0] == 0.001}, {x, y}, {t, 0, NL}]

Is this question the same as your previous post?

Get different Lyapunov exponent of DDE? Get different Lyapunov exponent of DDE?

No, my objective is to calculate Lyapunov exponent (LE) for delay differential equation(DDE) but the procedure is different in the two questions. There are lots of example for finding out LE of ordinary differential equation in mathematica but I am working with DDE and I am not getting anything regarding this. I have given my equation, please help me in this connection.