arXiv:1011.1508v1 [] 5 Nov 2010

Forecast Bias Correction: A Second-Order Method
Sean Crowell
Department of Mathematics
University of Oklahoma
Norman, OK 73019
scrowell@ou.edu

S. Lakshmivarahan*
School of Computer Science
University of Oklahoma
Norman, OK 73019
varahan@ou.edu

Abstract
In this paper we describe a methodology for determining corrections to
parameters and initial conditions in order to improve model forecasts in the
presence of data. This methodology is grounded in the variational problem
of minimizing the norm of the errors between forecast and data. The general
method is presented, and then the method is applied to the specific example
of a scalar model, the logistic equation, where it is shown that the method
produces satisfactory results.

Acknowledgement: We wish to record our sincere thanks to John Lewis for his interest in this work.

February 19, 2018

*Corresponding Author: varahan@ou.edu

1

Introduction

Models are abstractions of reality and the goodness of a model is often judged by its
ability to explain the observations which are the reflections of the underlying reality.
The difference between the actual observation and the model counterpart of the observation is known as the prediction error or forecast bias. This error/bias is a result
of either an ”inadequate” model or a result of misspecifications of initial/boundary
conditions and/or parameters in the model. Recently, Lakshmivarahan and Lewis
(2008) (LL(2008), hereafter) have developed a framework for estimating the error
in initial/boundary conditions and/or parameters that will account for the observed
forecast bias. Using the (first-order) sensitivity of the solution with respect to the
initial conditions and parameters and the observed forecast errors, this framework
recasts the bias estimation problem as an inverse problem, in particular as a linear least-squares problem. In LL(2008), they demonstrated the power of this idea
using a simplified model consisting of a system of three coupled nonlinear ordinary
differential equations that describes the evolution of the mixed-layer over the Gulf
of Mexico and actual observations obtained from an earlier field experiment. For
details, refer to LL(2008) and references therein.
In this paper we further analyze the power of this framework using a scalar model
for population growth known as the logistic model in continuous time. But instead of
the first-order method used in LL(2008), in this paper we compare the performance
of first-order and second-order methods. For simplicity it is assumed that the model
is adequate and the forecast errors are due only to errors in the initial condition
and the parameter. As will become evident, the first-order method leads to a linear
1

least-squares problem and the second-order method leads to a nonlinear least-squares
problem. Analysis of model inadequacies will be explored in a companion paper.
Section 2 contains a summary of the framework and is adapted from LL(2008).
Properties of the logistic model are described in Section 3. Numerical experiments
related to the estimation of errors using the first-order and second-order methods
along with a comparison of these methods are contained in Section 4. Concluding
observations are contained in Section 5.

2

The Framework - a Summary

Let x(t) ∈ R denote the state of a dynamical system at time t ∈ R+ where x(0) is
called the initial state. Let f : RxRxR+ → R and the state evolve according to a
scalar nonlinear ordinary differential equation given by

ẋ = f (x, α, t)

(2.1)

where α ∈ R is called the parameter. It is tacitly assumed that the function f in (2.1)
satisfies all of the conditions required for the existence, uniqueness, and smoothness
of the solution x(t) = x(t, α, x(0)) for all t ∈ R+ .
Let h : R → R be a smooth function that defines the model counterpart of the
observation
ztM = h(x(t))

2

(2.2)

Let zt be the actual observation obtained from the field measurement. Then
et = zt − ztM

(2.3)

denotes the prediction error or the forecast bias. It is assumed that the model
in (2.1) does not have any deficiencies and the forecast bias is mainly due to the
misspecification of the initial condition, x(0), and/or the parameter α.
Let ∆x be the actual change in x(t) induced by the perturbations δx(0) in x(0)
and δα in α. The goal is to find β = (δx(0), δα)T, the vector of perturbations such
that
zt − h(xt + ∆x) = 0.

(2.4)

Let ∆h be the actual change in h(x(t)) induced by the change ∆x in x(t). Then

h(x(t) + ∆x) = h(x(t)) + ∆h

(2.5)

Combining (2.3) - (2.5), our goal is to find β such that

et = ∆h

(2.6)

Recall that
∆h =

∞
X

δj h

(2.7)

j=1

where δ j h is called the j t h variation of h, the fraction of the induced change that is

3

attributable to the j th derivative of h and the change in ∆x in x(t). Similarly

∆x =

∞
X

δj x

(2.8)

j=1

where δ j x, called the j th variation of x, is the fraction of the total change in x(t)
that is attributable to the j th partial derivatives of x(t) with respect to α and x(0)
and β.
In practice we approximate the infinite sum by taking only the first k terms,
resulting in a k th -order approximation given by

∆h =

k
X

j

δ h and ∆x =

δ j x.

(2.9)

j=1

j=1

2.1

k
X

First-order Approximation

Setting k = 1, we get
∆h = δh and ∆x = δx.
Let Da (g) =

∂g
.
∂a

(2.10)

From first principles and (2.10), it follows that

δh = Dx (h)∆x = Dx (h)δx

(2.11)

δx = Dx(0) (x)δx(0) + Dα (x)δα.

(2.12)

and

4

The first derivative Da (g) is called the first-order sensitivity of g with respect to a
(Cacuci (2003) and Cruz(1973)). Setting
H1 = Dx(0) (x), H2 = Dα (x), H = [H1 , H2 ] ∈ R1x2 ,

(2.12) can be rewritten as
δx = Hβ.

(2.13)

Combining (2.13) and (2.11) with (2.6) we obtain an underdetermined linear leastsquares problem
Hβ = e

(2.14)

H = Dx (h)H ∈ R1x2

(2.15)

where

Suppose that there are N observations zt1 , zt2 , ... ztN available at times 0 ≤ t1 <
t2 < ... < tN . Then, at each time ti , 1 ≤ i ≤ N, we have the forecast error

eti = zti − h(xti )

(2.16)

H1 (ti ) = Dx(0) (x(ti )), H2 (ti ) = Dα (x(ti ))

(2.17)

Define

H(ti ) = [H1 (ti ), H2 (ti )] ∈ R1x2

5

and a diagonal matrix

D(h) = Diag{D1 (h), D2 (h), ..., DN (h)} ∈ RN xN

where
Di (h) = Dx (h(x(ti ))).
for simplicity in notation. Let


H1 (t1 )

H2 (t1 )







 H1 (t1 ) H2 (t1 ) 


HN =  .
 ∈ RN x2
.
 .

.
. 
 .


H1 (tN ) H2 (tN )

(2.18)

HN = D(h)HN

(2.19)

EN = (et1 , et2 , ..., etN )T ∈ RN

(2.20)

and

For this case of N observations, in place of (2.14), we obtain a linear least-squares
problem
HN β = EN

(2.21)

where HN ∈ RN x2 , β ∈ R2 , and EN ∈ RN .
The unknown β is obtained by minimizing
g1 (β) = kHN β − EN k2
6

(2.22)

where it is tacitly assumed that the matrix HN is of full rank. Under this assumption,
β is given by (Lewis, et al (2006))

−1 T


HNT HN
HN EN if N > 2



β=
HN−1EN
if N = 2





H T HN H T −1 EN if N < 2
N
N

(2.23)

When HN is not of full rank, we invoke the Tikhonov regularization (Lewis, et al
(2006)) using which β is obtained by minimizing
g1 (β) = kHN β − EN k2 +

λ2
kβk2
2

(2.24)

for some real constant λ > 0, called the regularization parameter. The minimizing
β in this case is given by

β LS = (HNT HN + λI)−1 HNT EN

(2.25)

In place of Tikhonov regularization, one could use the generalized inverse of the
matrix HN to obtain
β = HN+ EN
where HN+ is the Moore-Penrose inverse (Lewis et al (2006)).

7

2.2

Second-Order Approximation

Setting k = 2, in (2.9), it follows that
∆x = δx + δ 2 x

(2.26)

where δx is given by (2.13) is linear in β. From first principles, we obtain
1
δ 2 x = β T D 2 (x)β
2

(2.27)

which is quadratic in β and D 2 (x) is the Hessian of x(t) with respect to the entries
of β, that is




2
2
 Dx(0) (x(t)) Dx(0),α (x(t))
D 2 (x) = 

2
Dα,x(0)
(x(t))
Dα2 (x(t))

where Da2 (g) =

∂2g
∂a2

2
and Da,b
(g) =

∂2g
.
∂a∂b

(2.28)

The second partial derivatives of g are also

called the second-order sensitivities of g with respect to a and b. Similarly, from
(2.9),
1
∆h = δh + δ 2 h = Dx (h)∆x + Dx2 (h)(∆x)2
2

(2.29)

Substituting (2.26) in (2.29) and using (2.6) for the case of the single observation,
we obtain a nonlinear least-squares problem.
1
et = Dx (h)(δx + δ 2 x) + Dx2 (h)(δx + δ 2 x)2
2

8

(2.30)

The unknown β is then obtained by minimizing
1
g2 (β) = et − Dx (h)(δx + δ x) + Dx2 (h)(δx + δ 2 x)2
2

2

2

(2.31)

In the special case where the state is directly observable, h(x) = x and so Dx (h) = 1
and Dx2 (x) = 0. Substituting these into (2.31), in view of (2.13) and (2.17), we get
1
g2 (β) = et − Hβ − β T D 2 (x)β
2

2

(2.32)

Clearly, g2 (β) is a fourth degree polynomial in the components of β and in general
can have multiple minima which could further complicate the minimization problem.
When there are N (≥ 2) observations at times 0 ≤ t0 < t1 < · · · < tN −1 , then we get
N versions of the relation (2.30), one for each time. The unknown β is then obtained
by minimizing

G2 (β) =

N
−1
X
i=0

1 2
(h)[δx(ti ) + δ 2 x(ti )]2
eti − Dx(ti ) (h)[δx(ti ) + δ x(ti )] − Dx(t
i)
2
2

2

.
(2.33)

In the special case when h(x) = x, Dx (h) = 1 and Dx2 (x) = 0 and (2.33) reduces to

G2 (β) =

N
−1
X
i=0

1
eti − H(ti )β − β T D 2 (x(ti )β
2

9

2

.

(2.34)

3

The Logistic Equation

Consider the standard logistic equation that describes the growth of a certain population in an environment with carrying capacity 1 and growth rate parameter α > 0
given by
ẋ = αx(1 − x)

(3.1)

with x(0) = xo > 0 as the initial condition. The solution of this nonlinear ordinary
differential equation is given by

x(t) =

xo eαt
1 − xo + xo eαt

(3.2)

Notice that x(t) depends nonlinearly on xo and α.
It can be verified that x(t) → 1 as t → ∞ and x(t) ∈ [0, 1] for all t whenever
x(0) ∈ [0, 1], for α > 0. Since ẋ(t) > 0 for all t, and so is an increasing function
which increases from xo to 1 as t increases from 0 to ∞. A typical plot of x(t) in
(3.2) for x(0) = 0.5 and α = 1 is given in Figure 1.
The various first-order sensitivities of the the solution with respect to β =
(x(0), α)T are given by
eαt
[1 − xo + xo eαt ]2
xo (1 − xo )teαt
Dα (x) =
= xo (1 − xo )tDx(0) (x)
[1 − xo + xo eαt ]2

Dx(0) (x) =

10

(3.3)
(3.4)

That is, Dα (x) is a multiple of Dx(0) . Further it can be verified for large t that

Dα (x) ≈

1
1 − xo −αt
te
and Dx(0) (x) ≈ 2 e−αt
xo
xo

Thus, Dα (x) and hence Dx(0) (x) tend to zero as t → ∞. A plot of Dα (t) and
Dx(0) (x) versus t for xo = 0.5 and α = 1 is given in Figure 2. This in turn implies
that x(t) becomes less sensitive with respect to α and x(0) for large t, and that for
our framework to be effective, the observations have to be taken during the transient
phase of the solution.
The Hessian of x(t) with respect to β whose elements are the second-order sensitivity functions are given by




2
2
 Dx(0) (x) Dx(0),α (x)
D 2 (x) = 

2
Dα,x(0)
(x)
Dα2 (x)

(3.5)

where
2
(x) =
Dx(0)

Dα2 (x) =

2(eαt − e2αt )
[1 − xo + xo eαt ]3

xo (1 − xo )t2 eαt [1 − xo − xo eαt ]
[1 − xo + xo eαt ]3

and
2
2
Dx(0),α
(x) = Dα,x(0)
(x) =

t2 eαt [1 − xo − xo eαt ]
.
[1 − xo + xo eαt ]3

It is assumed that x(t) is directly observable, and so h(x) = x. Accordingly Dx (h) ≡
1 and Dx2 (h) ≡ 0.

11

3.1

The First-Order Approximation

In this case (2.14) reduces to (2.13) and from (3.3) and (3.4) we get
H = H = Dx(0) (x)[1, x0 (1 − x0 )t] ∈ R1x2 .

(3.6)

Two cases arise. The first case is when there are less observations than parameters.
The second is when the number of observations is at least as great as the number
of unknown parameters. In this paper, we only address the latter case. In this case,
N > 1 and





D (x(t1 )) Dx(0) (x(t1 ))xo (1 − xo )t1

 x(0)


 Dx(0) (x(t2 )) Dx(0) (x(t2 ))xo (1 − xo )t2 


HN = 

.
.


..
..




Dx(0) (x(tN )) Dx(0) (x(tN ))xo (1 − xo )tN

(3.7)

It can be verified that matrix HN ∈ RN x2 is of rank two and the least-square solution
is given by
f
βLS

3.2

=





HN−1 EN

if N = 2

(3.8)


(HNT HN )−1 HNT EN if N > 2

The Second-Order Approximation

Expanding the right hand side of (2.34) and dropping the terms of degree 3 and 4,
we obtain a quadratic approximation Q(β) to g2 (β) given by

Q(β) =

N
−1
X
i=0


 2
eti − 2eti H(ti )β + β T [H(ti )T H(ti ) − eti D 2 (x(ti ))]β .

12

(3.9)

The β that minimizes this is given by

s
βLS

=

N
−1
X
i=0

4




H(ti )T H(ti ) − eti D 2 (x(ti )) H(ti )T eti .

(3.10)

Numerical Experiments

In this section we compare the performance of the first-order and the second-order
methods using N (≥ 2) observations.

4.1

Generation of Observations

It is assumed that h(x) = x, that is, the state of the system x(t) is directly observable
and that zt = x(t). We use the model solution x(t) starting from x(0) and α to
generate different sets of N observations.
The distribution of the N observation times

t0 < t1 < t2 < · · · < tN −1

are given in Table 1 where ti = t0 + i(k∆) with ∆ = 0.5 (fixed), k = 1, 4, 8 and
12 and 0 ≤ i ≤ N − 1. Here t0 denotes the starting time and (k∆) denotes the
time between successive observations. Thus, for example, for t0 = 0 and k = 4, we
generate N = 6 observations at times t0 = 0.0, t1 = 2.0, t2 = 4.0, t3 = 6.0, t4 = 8.0,
and t5 = 10.0, which corresponds to the second row of Table 1.
The actual observations zti = x(ti ) where x(ti ) is the solution of (3.1) starting
from x(0) = 0.5 and α = 1.0 used in the estimation of β are given in Table 2.
13

4.2

A Plan for the Experiment

Let x(t) be the solution of the model (3.1) starting from x(0) 6= x(0) and α 6= α. Our
goal is to estimate δx(0) = x(0) − x(0) and δα = α − α using the set of observations
in Table 2. We divide the experiment into three parts based on N, the number of
observations.
Case 1: N = 2. We compute the estimates of β using two observations at times t0
and t0 + k∆ for k = 1, 4, 8 and 12 and for t0 = 0, 4 and 8. Estimates of β using the
first-order and the second-order methods are given in the top part of Tables 3 and 4
respectively.
For example, referring to the top part of Table 3, the estimate of β = (δx(0) δα)
using two observations zt0 = z4 = x(4) and zt1 = z4.5 = x(4.5) and the first-order
method is given by β = (−0.0955, 0.0981). Since the actual value of δx(0) = −0.1
and δα = 0.1, the relative error in the estimate of δx(0) =
and that of δα =

0.0981−0.1
0.1

−0.0955−(−0.1)
−0.1

= 0.045

= −0.019. Similarly for all other patterns of N = 2

observations.
Case 2: N = 4. Results of the estimation of β using 4 observations at times t0 ,
t0 + k∆, t0 + 2k∆ and t0 + 3k∆ for k = 1, 4, 8 and 12 and t0 = 0, 4 and 8 are given
in the middle parts of Tables 3 and 4.
Case 3: N = 6. Results of the estimation of β using 4 observations at times t0 ,
t0 + k∆, t0 + 2k∆, ..., t0 + 5k∆ for k = 1, 4, 8 and 12 and t0 = 0, 4 and 8 are given
in the bottom parts of Tables 3 and 4.

14

4.3

Comments and Discussion

A number of observations are in order.
1. When two observations are too close to each other, say by a distance ∆ = 0.5,
the solutions at these two times are very close to each other and consequently
the sensitivity functions at time t0 (Dx(0) (x(t0 )), Dα (x(t0 ))) and at time t1 =
t0 + ∆ namely (Dx(0) (x(t1 )), Dα (x(t1 ))) are also very close to each other. This
in turn implies that the two rows in H N are very nearly colinear, which leads
T

to bad conditioning of the matrix H N H N . Hence there is a greater chance for
the estimates to be less reliable.
2. From (3.2) and Figure 1, it follows that the solution x(t) attains the steady state
value which is independent of x(0) and α. In other words, the sensitivity of x(t)
with respect to x(0) and α decrease to zero as t increases. Thus, from (3.3) and
(3.4) it follows that Dx(0) (x(t)) and Dα (x(t)) both tend to zero as t increases.
Thus, the row of H N for large t becomes zero leading to ill-conditioning of H N
which makes the estimates less reliable.
3. Since the solution x(t) depends on x(0) and α nonlinearly, we cannot hope to
recover the actual perturbations δx0 and δα that led to the forecast bias in
the first place. While in principle the second-order method is better than the
first-order counterpart, to improve the overall accuracy of the estimate, one
may have to resort to an iterative improvement of the estimate. A detailed
comparison of the performance of the iterative versions of the first and second
order methods for small values of the perturbations (δx0 = −0.1 and δα = 0.1)
15

are given in Tables 5 and 6 respectively. Tables 7 and 8 contain similar results
for large values of δx0 and δα.
It turns out that the iterative versions of the first-order method peform at least
as well as the iterative versions of the second-order counterpart. Given that
second-order methods require larger computation, from this exercise it follows
that the iterative first-order method would be a good choice for estimating the
forecast bias.

16

5

References

Cacuci, D.G. (2003) Sensitivity and Uncertainty Analysis, Chapman and Hall/CRC
Press, Boca Raton, FL

Cruz, Jr. J.B. (1973) System Sensitivity Analysis, Stroudsberg, PA. (editor)

Lakshmivarahan, S. and J. Lewis (2008) ”Finding sources of bias error in forecast
models: A framework,” Monthly Weather Review. (submitted)

Lewis, J., S. Lakshmivarahan, and S.K. Dhall (2006) Dynamic Data Assimilation,
Cambridge University Press, 654 pages

17

18

Figure 1: Typical Solution of Logistic Equation (x(0) = 0.5 and α = 1)

(a) Sensitivity with respect to x(0)

(b) Sensitivity with respect to α

Figure 2: First-Order Sensitivities of the Solution to the Logistic Equation

19

k∆
t0
t1
t2
t3
t4
t5
k∆
t0
t1
t2
t3
t4
t5
k∆
t0
t1
t2
t3
t4
t5

Observation
∆
4∆
0.0 0.0
0.5 2.0
1.0 4.0
1.5 6.0
2.0 8.0
2.5 10.0
∆
4∆
4.0 4.0
4.5 6.0
5.0 8.0
5.5 10.0
6.0 12.0
6.5 14.0
∆
4∆
8.0 8.0
8.5 10.0
9.0 12.0
9.5 14.0
10.0 16.0
10.5 18.0

Times
8∆
0.0
4.0
8.0
12.0
16.0
20.0
8∆
4.0
8.0
12.0
16.0
20.0
24.0
8∆
8.0
12.0
16.0
20.0
24.0
28.0

12∆
0.0
6.0
12.0
18.0
24.0
30.0
12∆
4.0
10.0
16.0
22.0
28.0
34.0
12∆
8.0
14.0
20.0
26.0
32.0
38.0

Table 1: Distribution of the N observation times: ti = t0 + k(i − 1)∆ with ∆ = 0.5,
k = 1, 4, 8, and 12 and 0 ≤ i ≤ N − 1. t0 denotes the starting time and k∆ denotes
the time interval between successive observations and N denotes the number of
observations.

20

t0 = 0
z0
z1
z2
z3
z4
z5

1
0.50000000000000
0.62245933120185
0.73105857863000
0.81757447619364
0.88079707797788
0.92414181997876

t0 = 4
z0
z1
z2
z3
z4
z5

1
0.98201379003791
0.98901305736941
0.99330714907572
0.99592986228410
0.99752737684337
0.99849881774326

t0 = 8
z0
z1
z2
z3
z4
z5

1
0.99966464986953
0.99979657302194
0.99987660542401
0.99992515377249
0.99995460213130
0.99997246430889

Actual Observations
Multiple of ∆
4
8
0.50000000000000 0.50000000000000
0.88079707797788 0.98201379003791
0.98201379003791 0.99966464986953
0.99752737684337 0.99999385582540
0.99966464986953 0.99999988746484
0.99995460213130 0.99999999793885
Multiple of ∆
4
8
0.98201379003791 0.98201379003791
0.99752737684337 0.99966464986953
0.99966464986953 0.99999385582540
0.99995460213130 0.99999988746484
0.99999385582540 0.99999999793885
0.99999916847197 0.99999999996225
Multiple of ∆
4
8
0.99966464986953 0.99966464986953
0.99995460213130 0.99999385582540
0.99999385582540 0.99999988746484
0.99999916847197 0.99999999793885
0.99999988746484 0.99999999996225
0.99999998477002 0.99999999999931

Table 2: Actual Observations

12
0.50000000000000
0.99752737684337
0.99999385582540
0.99999998477002
0.99999999996225
0.99999999999991
12
0.98201379003791
0.99995460213130
0.99999988746484
0.99999999972105
0.99999999999931
1.00000000000000
12
0.99966464986953
0.99999916847197
0.99999999793885
0.99999999999489
0.99999999999999
1.00000000000000

21

True Corrections
2 observations
to
0
δxo
δα
κ
4
δxo
δα
κ
8
δxo
δα
κ
4 observations
to
0
δxo
δα
κ
4
δxo
δα
κ
8
δxo
δα
κ
6 observations
to
0
δxo
δα
κ
4
δxo
δα
κ
8
δxo
δα
κ

δxo = −0.1
δα = 0.1
multiple of ∆
1
4
8
12
-0.1000 -0.1000 -0.1000 -0.1000
0.0761 0.0968 0.1028 0.0989
2.8e+2 4.1e+1 2.0e+2 3.1e+3
-0.0955 -0.0888 -0.0809 -0.0739
0.0981 0.0912 0.0829 0.0756
1.3e+3 6.0e+2 5.0e+3 8.1e+4
-0.0479 -0.0390 -0.0284 -0.0190
0.0657 0.0611 0.0555 0.0507
7.9e+3 3.8e+3 3.2e+4 5.2e+5
multiple of ∆
1
4
8
12
-0.1006 -0.1000 -0.1000 -0.1000
0.0892 0.0977 0.1028 0.0989
7.0e+1 3.5e+1 2.0e+2 3.1e+3
-0.1038 -0.1108 -0.1477 0.4009
0.0977 0.1039 0.1349 -0.0570
3.9e+2 9.4e+2 5.0e+4 1.5e+5
-0.0440 -0.0378 -0.0283 -0.0190
0.0637 0.0605 0.0555 0.0507
2.5e+3 3.4e+3 3.2e+4 5.2e+5
multiple of ∆
1
4
8
12
-0.1011 -0.1000 -0.1000 -0.1000
0.0951 0.0977 0.1028 0.0989
4.2e+1 3.5e+1 2.0e+2 3.1e+3
-0.0907 -0.0879 -0.0808 -0.0739
0.0932 0.0902 0.0828 0.0756
3.0e+2 5.4e+2 5.0e+3 8.1e+4
-0.0416 -0.0378 -0.0283 -0.0190
0.0625 0.0605 0.0555 0.0507
1.8e+3 3.4e+3 3.2e+4 5.2e+5

Table 3: Estimates of β using the first-order method. κ denotes the condition number
T
of the matrix H N H N , where N is the number of observations used in the estimation.

22
True Corrections
2 observations
to
0
δxo
δα
κ
4
δxo
δα
κ
8
δxo
δα
κ
4 observations
to
0
δxo
δα
κ
4
δxo
δα
κ
8
δxo
δα
κ
6 observations
to
0
δxo
δα
κ
4
δxo
δα
κ
8
δxo
δα
κ

δxo = −0.1
δα = 0.1
multiple of ∆
1
4
8
12
-0.1000 -0.1000 -0.1000 -0.1000
0.1117 0.0987 0.0945 0.0980
2.1e+2 2.5e+1 1.8e+2 7.1e+3
-0.1009 -0.1098 -0.1489 -0.1264
0.0953 0.1030 0.1359 0.1247
1.1e+3 9.6e+2 5.1e+4 2.5e+5
0.3763 0.3744 0.3737 0.3749
-0.0528 -0.0523 -0.0521 -0.0524
8.2e+3 4.7e+3 5.2e+4 1.1e+6
multiple of ∆
1
4
8
12
-0.0995 -0.1000 -0.1000 -0.1000
0.1025 0.0983 0.0945 0.0980
5.4e+1 2.3e+1 1.8e+2 7.1e+3
-0.1038 -0.1108 -0.1477 0.4009
0.0977 0.1039 0.1349 -0.0570
3.9e+2 9.4e+2 5.0e+4 1.5e+5
0.3754 0.3743 0.3737 0.3749
-0.0525 -0.0523 -0.0521 -0.0524
2.6e+3 4.3e+3 5.2e+4 1.1e+6
multiple of ∆
1
4
8
12
-0.0991 -0.1000 -0.1000 -0.1000
0.0991 0.0983 0.0945 0.0980
3.3e+1 2.3e+1 1.8e+2 7.1e+3
-0.1059 -0.1108 -0.1477 0.4009
0.0995 0.1039 0.1349 -0.0570
3.0e+2 9.4e+2 5.0e+4 1.5e+5
0.3750 0.3743 0.3737 0.3749
-0.0524 -0.0523 -0.0521 -0.0524
1.9e+3 4.3e+3 5.2e+4 1.1e+6

Table 4: Estimates of β using the second-order method. κ denotes the condition
T
number of the matrix H N H N , where N is the number of observations used in the
estimation.

23
True Corrections
2 observations
t0
0
δx0
δα
κ
iterations
4
δx0
δα
κ
iterations
8
δx0
δα
κ
iterations
4 observations
to
0
δx0
δα
κ
iterations
4
δx0
δα
κ
iterations
8
δx0
δα
κ
iterations
6 observations
to
0
δx0
δα
κ
iterations
δx0
δα
κ
iterations
8
δx0
δα
κ
iterations

δxo = −0.1
δα = 0.1
multiple of ∆
1
4
8
12
-0.1000 -0.1000 -0.1000 -0.1000
0.1000
0.1000
0.1000
0.1000
2.6e+2 3.2e+1 2.0e+2 4.6e+3
4
4
4
4
-0.1000 -0.1000 -0.1000 -0.1000
0.1000
0.1000
0.1000
0.1000
1.4e+3 8.7e+2 1.1e+4 2.7e+5
4
4
4
5
-0.1000 -0.1000 -0.1000 -0.1000
0.1000
0.1000
0.1000
0.1000
8.6e+3 5.7e+3 7.5e+4 1.8e+6
5
5
5
5
multiple of ∆
1
4
8
12
-0.1000 -0.1000 -0.1000 -0.1000
0.1000
0.1000
0.1000
0.1000
6.3e+1 2.9e+1 2.0e+2 4.6e+3
4
4
4
4
-0.1000 -0.1000 -0.1000 -0.1000
0.1000
0.1000
0.1000
0.1000
4.6e+2 8.1e+2 1.1e+4 2.7e+5
4
4
4
5
-0.1000 -0.1000 -0.1000 -0.1000
0.1000
0.1000
0.1000
0.1000
2.9e+3 5.3e+3 7.4e+4 1.8e+6
5
5
5
5
multiple of ∆
1
4
8
12
-0.1000 -0.1000 -0.1000 -0.1000
0.1000
0.1000
0.1000
0.1000
3.9e+1 2.9e+1 2.0e+2 4.6e+3
4
4
4
4
-0.1000 -0.1000 -0.1000 -0.1000
0.1000
0.1000
0.1000
0.1000
3.5e+2 8.1e+2 1.1e+4 2.7e+5
4
4
4
5
-0.1000 -0.1000 -0.1000 -0.1000
0.1000
0.1000
0.1000
0.1000
2.2e+3 5.3e+3 7.4e+4 1.8e+6
5
5
5
5

Table 5: Estimates of β using the iterated first-order method. κ denotes the condition number of
T

the matrix H N H N , where N is the number of observations used in the estimation. The number of
iterations is that required to reach a threshold of 10− 6 (at most 10)

24
True Corrections
2 observations
t0
0
δx0
δα
κ
iterations
4
δx0
δα
κ
iterations
8
δx0
δα
κ
iterations
4 observations
to
0
δx0
δα
κ
iterations
4
δx0
δα
κ
iterations
8
δx0
δα
κ
iterations
6 observations
to
0
δx0
δα
κ
iterations
4
δx0
δα
κ
iterations
8
δx0
δα
κ
iterations

δxo = −0.1
δα = 0.1
multiple of ∆
1
4
8
12
-0.1000 -0.1000 -0.1000 -0.1000
0.1000
0.1000 0.1000
0.1000
2.6e+2 3.2e+1 2.0e+2 4.6e+3
4
4
4
4
-0.1000 -0.1000 -0.1000 -0.1000
0.1000
0.1000 0.1000
0.1000
1.4e+3 8.7e+2 1.1e+4 2.7e+5
4
4
5
6
-0.1000 -0.1000 -0.1000
NaN
0.1000
0.1000 0.1000
NaN
8.6e+3 5.7e+3 7.5e+4 1.0e+9
6
5
5
10
multiple of ∆
1
4
8
12
-0.1000 -0.1000 -0.1000 -0.1000
0.1000
0.1000 0.1000
0.1000
6.3e+1 2.9e+1 2.0e+2 4.6e+3
4
4
4
4
-0.1000 -0.1000 -0.1000 -0.0868
0.1000
0.1000 0.1000
0.1023
4.6e+2 8.1e+2 1.1e+4 2.4e+5
4
4
5
10
-0.1000 -0.1000 -0.1000 -0.1000
0.1000
0.1000 0.1000
0.1000
2.9e+3 5.3e+3 7.4e+4 1.8e+6
7
7
6
6
multiple of ∆
1
4
8
12
-0.1000 -0.1000 -0.1000 -0.1000
0.1000
0.1000 0.1000
0.1000
3.9e+1 2.9e+1 2.0e+2 4.6e+3
3
4
4
4
-0.1000 -0.1000 -0.1000
NaN
0.1000
0.1000 0.1000
NaN
3.5e+2 8.1e+2 1.1e+4 1.0e+12
4
4
5
10
-0.1000 -0.1000 -0.1000
NaN
0.1000
0.1000 0.1000
NaN
2.2e+3 5.3e+3 7.4e+4 7.2e+12
7
7
6
10

Table 6: Estimates of β using the iterated second-order method. κ denotes the condition number
T

of the matrix H N H N , where N is the number of observations used in the estimation. The number
of iterations is that required to reach a threshold of 10− 6 (at most 10)

25

t0
0

4

8

True Parameter
Model Values
x0
α
0.3
0.8
0.3
0.9
0.3
1.1
0.3
1.2
0.4
0.8
0.4
1.1
0.4
0.9
0.4
1.2
0.6
0.8
0.6
0.9
0.6
1.1
0.6
1.2
0.7
0.8
0.7
0.9
0.7
1.1
0.7
1.2
0.3
0.8
0.3
0.9
0.3
1.1
0.3
1.2
0.4
1.1
0.4
1.2
0.4
0.8
0.4
0.9
0.6
0.8
0.6
0.9
0.6
1.1
0.6
1.2
0.7
0.8
0.7
0.9
0.7
1.1
0.7
1.2
0.3
0.8
0.3
1.1
0.3
0.9
0.3
1.2
0.4
0.8
0.4
0.9
0.4
1.1
0.4
1.2
0.6
0.8
0.6
1.1
0.6
1.2
0.6
0.9
0.7
1.1
0.7
1.2
0.7
0.8
0.7
0.9

Values: xo = 0.5, α = 1.0
Corrections
δx0
δα
Iterations
0.2000
0.2000
5
0.2000
0.1000
5
0.2000 -0.1000
5
0.2000 -0.2000
5
0.1000
0.2000
5
0.1000 -0.1000
4
0.1000
0.1000
4
0.1000 -0.2000
4
-0.1000 0.2000
4
-0.1000 0.1000
4
-0.1000 -0.1000
5
-0.1000 -0.2000
5
-0.2000 0.2000
4
-0.2000 0.1000
5
-0.2000 -0.1000
5
-0.2000 -0.2000
5
0.2000
0.2000
7
0.2000
0.1000
6
0.2000 -0.1000
5
0.2000 -0.2000
5
0.1000 -0.1000
4
0.1000 -0.2000
5
0.1000
0.2000
6
0.1000
0.1000
5
-0.1000 0.2000
5
-0.1000 0.1000
4
-0.1000 -0.1000
6
-0.1000 -0.2000
8
-0.2000 0.2000
5
-0.2000 0.1000
5
-0.2000 -0.1000
7
-0.2000 -0.2000
10
0.2000
0.2000
7
0.2000 -0.1000
5
0.2000
0.1000
7
0.2000 -0.2000
7
0.1000
0.2000
7
0.1000
0.1000
6
0.1000 -0.1000
5
0.1000 -0.2000
11
-0.1000 0.2000
6
-0.1000 -0.1000
6
NaN
NaN
100
-0.1000 0.1000
5
-0.2000 -0.1000
9
NaN
NaN
100
-0.2000 0.2000
6
-0.2000 0.1000
5

Table 7: Iterated First-Order Corrections Using Different Model Parameters α and
Initial Conditions x0 , in all cases using 4 observations spaced 0.5 seconds apart.

26

t0
0

4

8

True Parameter
Model Values
x0
α
0.3
0.8
0.3
0.9
0.3
1.1
0.3
1.2
0.4
0.8
0.4
0.9
0.4
1.1
0.4
1.2
0.6
0.8
0.6
0.9
0.6
1.1
0.6
1.2
0.7
0.8
0.7
0.9
0.7
1.1
0.7
1.2
0.3
0.8
0.3
0.9
0.3
1.1
0.3
1.2
0.4
0.8
0.4
0.9
0.4
1.1
0.4
1.2
0.6
0.8
0.6
0.9
0.6
1.1
0.6
1.2
0.7
0.8
0.7
0.9
0.7
1.1
0.7
1.2
0.3
0.8
0.3
1.1
0.3
1.2
0.3
0.9
0.4
0.8
0.4
0.9
0.4
1.1
0.4
1.2
0.6
0.8
0.6
0.9
0.6
1.1
0.6
1.2
0.7
0.8
0.7
0.9
0.7
1.1
0.7
1.2

Values: xo = 0.5, α = 1.0
Corrections
δx0
δα
Iterations
0.2000
0.2000
5
0.2000
0.1000
5
0.2000 -0.1000
5
0.2000 -0.2000
5
0.1000
0.2000
5
0.1000
0.1000
5
0.1000 -0.1000
4
0.1000 -0.2000
4
-0.1000 0.2000
4
-0.1000 0.1000
4
-0.1000 -0.1000
4
-0.1000 -0.2000
4
-0.2000 0.2000
4
-0.2000 0.1000
4
-0.2000 -0.1000
5
-0.2000 -0.2000
5
NaN
NaN
101
0.2000
0.1000
5
0.2000 -0.1000
5
0.2000 -0.2000
7
0.1000
0.2000
7
NaN
NaN
100
0.1000 -0.1000
4
0.1000 -0.2000
5
-0.1000 0.2000
7
-0.1000 0.1000
4
-0.1000 -0.1000
5
-0.1000 -0.2000
5
-0.2000 0.2000
5
-0.2000 0.1000
4
-0.2000 -0.1000
5
-0.2000 -0.2000
5
NaN
NaN
100
0.2000 -0.1000
6
0.2000 -0.2000
5
0.2000
0.1000
16
NaN
NaN
100
NaN
NaN
100
0.1000 -0.1000
5
0.1000 -0.2000
4
NaN
NaN
100
-0.1000 0.1000
7
-0.1000 -0.1000
5
-0.1000 -0.2000
6
-0.2000 0.2000
12
-0.2000 0.1000
5
-0.2000 -0.1000
5
-0.2000 -0.2000
9

Table 8: Iterated Second-Order Corrections Using Different Model Parameters α and
Initial Conditions x0 , in all cases using 4 observations spaced 0.5 seconds apart.

