**What is Data 2 Linear Dynamics (Data2LD)?**

Data 2 Linear Dynamics (Data2LD) estimates the solution and the parameters of linear dynamical systems from incomplete and noisy observations of the underlying processes.

It uses a linear combination of spline basis functions to approximate the implicitly defined solution of the dynamical system. While also estimating the systems' parameters by requiring that this approximating solution adheres to the data.

**Papers:**

Carey, M., Ramsay, J.O. (2020) 'Fast Stable Parameter Estimation for Linear Dynamical Systems'. Computational Statistics & Data Analysis

Carey, M., Gath, E., Hayes, K. (2016) 'A generalized smoother for linear ordinary differential equations'. *Journal of Computational and Graphical Statistics*.

Carey, M., Gath, E., Hayes, K. (2014) 'Frontiers in financial dynamics'. Research in International Business and Finance

Let's consider a study on traumatic brain injury (TBI), which contributes to just under a third (30.5%) of all injury-related deaths in the US and is caused by a blow to the head. Figure (1) shows the 133 accelerometer readings taken over 55.2 milliseconds. The dashed line represents the impulse function which denotes the blow to the head.

**Figure (1)**

The laws of motion tell us that the acceleration x(t) can be modelled by a second-order linear ordinary differential equation (ODE) with input a unit pulse u(t) representing the blow to the head and shown in the dashed lines in Figure (1).

This ODE

\begin{equation}

\frac{\textrm{d}^2x(t)}{\textrm{d}t^2} + \beta_{0} x(t) + \beta_{1} \frac{\textrm{d}x(t)}{\textrm{d}t} + \alpha u(t)

\end{equation}

contains three parameters $\beta_{0},\beta_{1}$ and $\alpha,$ and these convey the rate of the restoring force (as $t \rightarrow \infty,$ the acceleration will tend to revert back zero), the rate of the friction force (as $t \rightarrow \infty,$ the oscillations in the acceleration reduce to zero) and the rate of force from the unit pulse.

While there are several methods for estimating ODE parameters with partially observed data, they are invariably subject to several problems including high computational cost, sensitivity to initial values or large sampling variability.

**D ata to linear dynamics** (

*) overcomes these issues and produces estimates of the ODE parameters that have less bias, a smaller sampling variance and a ten-fold improvement in computation in comparsion to existing approaches.*

**Data2LD**# Demonstration of Data2LD: Head Impact Analysis

This page provides a detailed description of the MATLAB calculations necessary to run Data2LD code for estimating the parameters and the solution of linear differential equations.

## Contents

- Set-up
- Load the data
- The second order linear differential equation
- Set up basis functions for X
- Create a functional data object for the external/forcing function
- Run Data2LD
- Plot the fitted curve, the confidence and prediction intervals
- Plot the estimated parameters of the ODE and their approximated CI with respect to rho
- Check the residuals
- Check the normality
- Check the fitted vrs residuals plot

## Set-up

```
clear
clear functions
```

## Load the data

load motorcycledata.txt motot = motorcycledata(:,2); % time in milliseconds motoy = motorcycledata(:,3); % deformation in 0.1 millimeters % Adjust the data for baseline and plot impact = 14; % impact time baseind = motot < impact; basey = mean(motoy(baseind)); % Remove baseline, change time, and convert to centimeters motoy = (basey - motoy)/100.0; n = length(motoy); % Create a cell containing the data yCell = cell(1); yCell{1}(:,1) = motot; yCell{1}(:,2) = motoy; % Plot the data figure(1) phdl=plot(yCell{1}(:,1), yCell{1}(:,2), 'ko', [0,60], [0,0], 'k:'); set(phdl, 'LineWidth', 2) axis([min(yCell{1}(:,1))-0.01,max(yCell{1}(:,1))+0.01,-1.0,1.5]) xlabel('fontsize{16} Time (milliseconds)') ylabel('fontsize{16} Acceleration (cm/msec^2)')

Figure (1) illustrates observations of head acceleration measured 14 milliseconds before and 42.6 milliseconds after a blow to the cranium of a cadaver. The dashed line represents the impulse function which denotes the blow to the cranium that lasted 1 millisecond. The experiment, a simulated motor-cycle crash, is described in detail by Schmidt et al (1981).

## The second order linear differential equation

The second order linear differential equation, as used in this example is

\begin{equation}

\frac{\textrm{d}^2x(t)}{\textrm{d}t^2} + \beta_{0} x(t) + \beta_{1} \frac{\textrm{d}x(t)}{\textrm{d}t} + \alpha u(t)

\end{equation}

The three parameters , in the linear ODE convey the period of the oscillation, the change in its amplitude, as the oscillations decay exponentially to zero, and the size of the impact from the unit impulse respectively.

Our aim is to estimate the acceleration and the parameters from the data in Figure (1).

## Set up basis functions for X

motorng = [yCell{1}(1,1),yCell{1}(end,1)]; knots = [motorng(1),14,14,14,15,15,linspace(15,motorng(end),11)]; norder = 5; nbasis = length(knots)+ (norder - 2); Xbasisobj = create_bspline_basis(motorng,nbasis,norder,knots); %Create a cell for the basis for each variable XbasisCell = cell(1); XbasisCell{1} = Xbasisobj; %Plot the basis figure() plot(Xbasisobj)

For the expansion of we used order five B-spline functions, which by their nature have discontinuous third derivatives if all knots are singletons. To achieve curvature discontinuity at the impact point and at that point plus one, we placed three knots at these locations. We put no knots between the first observation and the impact point, where the data indicate a flat trajectory, and eleven equally spaced knots between the impact point plus one and the final time of observation.

## Create a functional data object for the external/forcing function

%The known external function is a unit impulse Ubasis = create_bspline_basis(motorng, 3, 1, [motorng(1),14,15,motorng(2)]); Ufd = fd([0;1;0],Ubasis); %Plot the external function figure() plot(Ufd) %Create a cell for the external functions EFCell = cell(1); EFCell{1} = Ufd;

Three order one B-splines over the knots [0, 14, 15, 56] with coefficient vector [0,1,0] represented the unit pulse function

## Run Data2LD

%Intial parameter estimates intial_theta = {0.1,0.1,0.1}; %Which parameters would we like to estimate. In this case all parameters. estimate = [1,1,1]; %Run Data2LD [Res,Results_cell] = D2LD_ODE(yCell,XbasisCell,intial_theta,... estimate,EFCell);

Elapsed time is 0.159864 seconds. Elapsed time is 0.158650 seconds. Elapsed time is 0.013674 seconds. Rho = 0.5 Iter. Criterion Grad Length 0 0.045955 0.001943 1 0.045726 0.000386 2 0.045676 0.000061 3 0.045673 0.000008 4 0.045673 0.000001 5 0.045673 0.000000 Convergence reached Rho = 0.73106 Iter. Criterion Grad Length 0 0.045718 0.000116 1 0.045716 0.000054 2 0.045715 0.000002 3 0.045715 0.000000 Convergence reached Rho = 0.8808 Iter. Criterion Grad Length 0 0.045854 0.000433 1 0.045851 0.000265 2 0.045841 0.000085 3 0.045841 0.000017 4 0.045841 0.000006 5 0.045841 0.000001 Convergence reached Rho = 0.95257 Iter. Criterion Grad Length 0 0.046300 0.002776 1 0.046286 0.000403 2 0.046281 0.000120 3 0.046281 0.000033 4 0.046281 0.000011 5 0.046281 0.000003 Convergence reached Rho = 0.97069 Iter. Criterion Grad Length 0 0.046799 0.003304 1 0.046789 0.000611 2 0.046780 0.000160 3 0.046780 0.000041 4 0.046779 0.000010 5 0.046779 0.000002 Convergence reached Rho = 0.98201 Iter. Criterion Grad Length 0 0.047602 0.006326 1 0.047583 0.001372 2 0.047556 0.000473 3 0.047553 0.000150 4 0.047553 0.000041 5 0.047553 0.000012 6 0.047553 0.000003 Convergence reached Rho = 0.98901 Iter. Criterion Grad Length 0 0.048725 0.010017 1 0.048698 0.002738 2 0.048626 0.000630 3 0.048623 0.000004 4 0.048623 0.000000 Convergence reached Rho = 0.99331 Iter. Criterion Grad Length 0 0.050119 0.013939 1 0.050085 0.005582 2 0.049909 0.000710 3 0.049906 0.000099 4 0.049906 0.000002 5 0.049906 0.000000 Convergence reached Rho = 0.99593 Iter. Criterion Grad Length 0 0.051423 0.017094 1 0.051388 0.010112 2 0.051152 0.002012 3 0.051150 0.000102 4 0.051150 0.000010 5 0.051150 0.000004 Convergence reached Rho = 0.99753 Iter. Criterion Grad Length 0 0.052353 0.016602 1 0.052329 0.013570 2 0.052175 0.002229 3 0.052175 0.000187 4 0.052175 0.000086 5 0.052175 0.000007 Convergence reached Rho = 0.9985 Iter. Criterion Grad Length 0 0.053037 0.013911 1 0.053024 0.015415 2 0.052946 0.001457 3 0.052946 0.000280 4 0.052946 0.000045 Convergence reached Rho = 0.99909 Iter. Criterion Grad Length 0 0.053532 0.010652 1 0.053525 0.014624 2 0.053490 0.000924 3 0.053490 0.000239 4 0.053490 0.000062 Convergence reached Rho = 0.99945 Iter. Criterion Grad Length 0 0.053873 0.007592 1 0.053870 0.010930 2 0.053855 0.000854 3 0.053855 0.000173 4 0.053855 0.000068 Convergence reached Rho = 0.99966 Iter. Criterion Grad Length 0 0.054099 0.005132 1 0.054097 0.006063 2 0.054091 0.001206 3 0.054091 0.000146 4 0.054091 0.000067 Convergence reached Rho = 0.9998 Iter. Criterion Grad Length 0 0.054243 0.003350 1 0.054243 0.003442 2 0.054241 0.001068 3 0.054241 0.000109 4 0.054241 0.000053 Convergence reached Rho = 0.99988 Iter. Criterion Grad Length 0 0.054335 0.002153 1 0.054335 0.001837 2 0.054334 0.000916 3 0.054334 0.000086 4 0.054334 0.000041 Convergence reached Rho = 0.99993 Iter. Criterion Grad Length 0 0.054392 0.001454 1 0.054392 8.390675 Rho = 0.99995 Iter. Criterion Grad Length 0 0.054430 0.002710 1 0.054430 12.076390 Elapsed time is 2.639805 seconds. Optimal values: Rho = 1 Theta1 = -0.057 Theta2 = -0.155 Theta3 = 0.402 sdTheta1 = 0.003 sdTheta2 = 0.019 sdTheta3 = 0.032 MSE = 0.054 ISE = 0 Df = 2.234 GCV = 0.056

## Plot the fitted curve, the confidence and prediction intervals

figure() hp = patch([yCell{1}(:,1); yCell{1}(end:-1:1,1);yCell{1}(1)],... [Results_cell{4}(:,4); Results_cell{4}(end:-1:1,5); Results_cell{4}(1)],[0.85 0.85 0.85]); hold on; plot(yCell{1}(:,1),[Results_cell{4}(:,2),Results_cell{4}(:,3)],'k--') plot(yCell{1}(:,1),Results_cell{4}(:,1),'k-') plot(yCell{1}(:,1), yCell{1}(:,2), 'ko'); xlim([min(yCell{1}(:,1))-0.1,max(yCell{1}(:,1))+0.1]) ylim([min(Results_cell{4}(:,4))-0.1,max(Results_cell{4}(:,5))+0.1]) xlabel('fontsize{16} Time (milliseconds)') ylabel('fontsize{16} Acceleration (cm/msec^2)') legend('Prediction Interval','Confidence Interval','Confidence Interval','Fitted Curve','Data') set(findall(gca, 'Type','text'), 'FontSize', 20)

This Figure illustrates the accelerometer readings of the brain tissue before and after a blow to the cranium are indicated by the circles. The fitted curve produced by Data2LD with $hat{rho} = 0.99$ (solid line), the approximated 95% pointwise confidence interval for the curve (dashed line) and the approximated 95% pointwise prediction interval for the curve (grey region)

## Plot the estimated parameters of the ODE and their approximated CI with respect to rho

Res_Mat = table2array(Res); for i=1:3 subplot(3,1,i) indexOfInterest = (Res_Mat(:,1) > 0.95); plot(Res_Mat(indexOfInterest,1),Res_Mat(indexOfInterest,i+1),'-*k') hold on; errorbar(Res_Mat(indexOfInterest,1),Res_Mat(indexOfInterest,i+1),... t_v.*Res_Mat(indexOfInterest,i+4),'ko') xlabel('$rho

The values of the three parameters and their approximated 95% confidence intervals for the head impact data over values of $rho$ converging to 1.0. Note the reduction in the width of the approximated 95% confidence intervals as $rho$ increase and the overlap of the final estimates of the ODE parameters, which illustrates that the parameter values have stabilised for high values of $rho$.

## Check if the residuals are approx. normal

figure() h = normplot(yCell{1}(:,2)-Results_cell{4}(:,1)); h(1).LineWidth = 2; h(2).LineWidth = 2; h(3).LineWidth = 2; % Our qqplot looks pretty good and indicates a normal distribution, as it’s generally a straight line.

## Check the fitted vrs residuals plot

figure() plot(Results_cell{4}(:,1),yCell{1}(:,2)-Results_cell{4}(:,1),'ko') % Should see random scatter in a constant band around zero that is no patterns, and no funnels.