Chapter 30. The Kalman Filter
270
matrix S = ksmooth()
# smoothed state only
matrix P
S = ksmooth(&P)
# the variance is wanted
These values are computed via a backward pass of the filter, from t  T to t 1, as follows:
L
t
F
t
K
t
H
0
t
u
t 1
H
t
Ö
Ö
Ö
1
t
e
t
L
0
t
u
t
U
t 1
H
t
ÖÖÖ
1
t
H
0
t
L
0
t
U
t
L
t
ˆ

tjT
ˆ

tjt 1
P
tjt 1
u
t 1
P
tjT
P
tjt 1
P
tjt 1
U
t 1
P
tjt 1
with initial values u
T
0 and U
T
0.
6
This iteration is preceded by a special forwardpass in which the matrices K
t
,ÖÖÖ
1
t
,
ˆ

tjt 1
and P
tjt 1
are stored for all t. If F is time-varying, its values for all t are stored on the forward pass, and
similarly for H.
30.8 The ksimul function
This simulation function takes up to three arguments. The first, mandatory, argument is a (T r)
matrix containing artificial disturbances for the state transition equation: row t of this matrix
represents v
0
t
. If the current filter has a non-null R (obsvar) matrix, then the second argument
should be a (T n) matrix containing artificial disturbances for the observation equation, on the
same pattern. Otherwise the second argument should be given as null. If r  1 you may give a
series for the first argument, and if n 1 a series is acceptable for the second argument.
Provided that the current filter does not include exogenous variables in the observation equation
(obsx), theT for simulation neednot equal that definedby the original obsydata matrix: in effect T
is temporarily redefinedby the rowdimension of thefirst argument to ksimul. Once the simulation
is completed, the T value associated with the original data is restored.
The value returned by ksimul is a (T n) matrix holding simulated values for the observables at
each time step. A third optional matrix-pointer argument allows you to retrieve a (T r) matrix
holding the simulated state vector. Examples:
matrix Y = ksimul(V)
# obsvar is null
Y = ksimul(V, W)
# obsvar is non-null
matrix S
Y = ksimul(V, null, &S) # the simulated state is wanted
The initial value 
1
is calculated thus: we find the matrix T such that TT0  P
1j0
(as given by the
inivar element in the kalman block), multiply it into v
1
, and add the result to 
1j0
(as given by
inistate).
If the disturbances are correlated across the two equations the arguments to ksimul must be
revised: the first argument should be a (T  p) matrix, each row of which represents "
"
"
0
t
(see sec-
tion30.2), and the second argument should be given as null.
30.9 Example 1: ARMA estimation
As iswell known, the Kalman filter providesa very efficient way to compute the likelihoodofARMA
models; as an example, take an ARMA(1,1) model
y
t
y
t 1
"
t
"
t 1
6
SeeI.Karibzhanov’sexpositionathttp://karibzhanov.com/help/kalcvs.htm.
Add pdf to powerpoint - SDK application project:C# Create PDF from PowerPoint Library to convert pptx, ppt to PDF in C#.net, ASP.NET MVC, WinForms, WPF
Online C# Tutorial for Creating PDF from Microsoft PowerPoint Presentation
www.rasteredge.com
Add pdf to powerpoint - SDK application project:VB.NET Create PDF from PowerPoint Library to convert pptx, ppt to PDF in vb.net, ASP.NET MVC, WinForms, WPF
VB.NET Tutorial for Export PDF file from Microsoft Office PowerPoint
www.rasteredge.com
Chapter 30. The Kalman Filter
271
One of the ways the above equation can be cast in state-space form is by defining a latent process
t
1  L
1
"
t
. The observation equation corresponding to (30.2) is then
y
t

t

t 1
(30.11)
and the state transition equation corresponding to (30.1) is
"
t
t 1
#
"
0
1
0
#"
t 1
t 2
#
"
"
t
0
#
The gretl syntax for a corresponding kalman block would be
matrix H = {1; theta}
matrix F = {phi, 0; 1, 0}
matrix Q = {s^2, 0; 0, 0}
kalman
obsy y
obsymat H
statemat F
statevar Q
end kalman
Note that the observation equation (30.11) does not include an “error term”; this is equivalent
to saying that Vw
t
  0 and, as a consequence, the kalman block does not include an obsvar
keyword.
Once the filter is set up, all it takes to compute the log-likelihood for given values of ,  and
2 is to execute the kfilter() function and use the $kalman_lnl accessor (which returns the
total log-likelihood) or, more appropriately if the likelihood has to be maximized through mle, the
$kalman_llt accessor, which returns the series of individual contribution to the log-likelihood for
each observation. An example is shown in script30.1.
30.10 Example 2: local level model
Suppose we have a series y
t
 
t
"
t
, where 
t
is a random walk with normal increments of
variance 
2
1
and "
t
is a normal white noise with variance 
2
2
,independent of 
t
. This is known as
the “local level” model in Harvey’s (1989) terminology, and it can be cast in state-space form as
equations (30.1)-(30.2) with F  1, v
t
N0; 
2
1
, H  1 and w
t
N0; 
2
2
. The translation to a
kalman block is
kalman
obsy y
obsymat 1
statemat 1
statevar s2
obsvar s1
end kalman --diffuse
The two unknown parameters 
2
1
and 
2
2
can be estimated via maximum likelihood. Script30.2
provides an example of simulation and estimation of such a model. For the sake of brevity, simu-
lation is carried out via ordinary gretl commands, rather than the state-space apparatus described
above.
The example contains two functions: the first one carries out the estimation of the unknown pa-
rameters 
2
1
and 
2
2
via maximum likelihood; the second one uses these estimates to compute
asmoothed estimate of the unobservable series 
t
under the name muhat. A plot of 
t
and its
estimate is presented in Figure30.1.
SDK application project:C# PDF insert image Library: insert images into PDF in C#.net, ASP
C#.NET PDF SDK - Add Image to PDF Page in C#.NET. How to Insert & Add Image, Picture or Logo on PDF Page Using C#.NET. Add Image to PDF Page Using C#.NET.
www.rasteredge.com
SDK application project:VB.NET PDF insert image library: insert images into PDF in vb.net
VB.NET PDF - Add Image to PDF Page in VB.NET. Guide VB.NET Programmers How to Add Images in PDF Document Using XDoc.PDF SDK for VB.NET.
www.rasteredge.com
Chapter 30. The Kalman Filter
272
Example 30.1: ARMA estimation
function void arma11_via_kalman(series y)
/* parameter initalization */
phi = 0
theta = 0
sigma = 1
/* Kalman filter setup */
matrix H = {1; theta}
matrix F = {phi, 0; 1, 0}
matrix Q = {sigma^2, 0; 0, 0}
kalman
obsy y
obsymat H
statemat F
statevar Q
end kalman
/* maximum likelihood estimation */
mle logl = ERR ? NA : $kalman_llt
H[2] = theta
F[1,1] = phi
Q[1,1] = sigma^2
ERR = kfilter()
params phi theta sigma
end mle -h
end function
# ------------------------ main ---------------------------
open arma.gdt
# open the "arma" example dataset
arma11_via_kalman(y) # estimate an arma(1,1) model
arma 1 1 ; y --nc
# check via native command
SDK application project:VB.NET PDF Password Library: add, remove, edit PDF file password
VB: Add Password to PDF with Permission Settings Applied. This VB.NET example shows how to add PDF file password with access permission setting.
www.rasteredge.com
SDK application project:C# PDF Password Library: add, remove, edit PDF file password in C#
C# Sample Code: Add Password to PDF with Permission Settings Applied in C#.NET. This example shows how to add PDF file password with access permission setting.
www.rasteredge.com
Chapter 30. The Kalman Filter
273
Example 30.2: Local level model
function matrix local_level (series y)
/* starting values */
scalar s1 = 1
scalar s2 = 1
/* Kalman filter set-up */
kalman
obsy y
obsymat 1
statemat 1
statevar s2
obsvar s1
end kalman --diffuse
/* ML estimation */
mle ll = ERR ? NA : $kalman_llt
ERR = kfilter()
params s1 s2
end mle
return s1 ~ s2
end function
function series loclev_sm (series y, scalar s1, scalar s2)
/* return the smoothed estimate of \mu_t */
kalman
obsy y
obsymat 1
statemat 1
statevar s2
obsvar s1
end kalman --diffuse
series ret = ksmooth()
return ret
end function
/* -------------------- main script -------------------- */
nulldata 200
set seed 202020
setobs 1 1 --special
true_s1 = 0.25
true_s2 = 0.5
v = normal() * sqrt(true_s1)
w = normal() * sqrt(true_s2)
mu = 2 + cum(w)
y = mu + v
matrix Vars = local_level(y)
# estimate the variances
muhat = loclev_sm(y, Vars[1], Vars[2]) # compute the smoothed state
SDK application project:C# PDF Sticky Note Library: add, delete, update PDF note in C#.net
C#.NET PDF SDK - Add Sticky Note to PDF Page in C#.NET. Able to add notes to PDF using C# source code in Visual Studio .NET framework.
www.rasteredge.com
SDK application project:C# WinForms Viewer: Load, View, Convert, Annotate and Edit
allowed to load and view PowerPoint without Microsoft Office software installed, create PDF file, Tiff image and HTML file from PowerPoint, add annotations to
www.rasteredge.com
Chapter 30. The Kalman Filter
274
-8
-6
-4
-2
0
2
4
6
8
10
0
50
100
150
200
mu
muhat
Figure 30.1: Local level model: 
t
and its smoothed estimate
By appending the following code snippet to the example in Table30.2, one may check the results
against the R commandStructTS.
foreign language=R --send-data
y <- gretldata[,"y"]
a <- StructTS(y, type="level")
a
StateFromR <- as.ts(tsSmooth(a))
gretl.export(StateFromR)
end foreign
append @dotdir/StateFromR.csv
ols muhat 0 StateFromR --simple
SDK application project:C# HTML5 Viewer: Load, View, Convert, Annotate and Edit PowerPoint
load and view PowerPoint without Microsoft Office software installed, convert PowerPoint to PDF file, Tiff image and HTML file, as well as add annotations in
www.rasteredge.com
SDK application project:RasterEdge XDoc.PowerPoint for .NET - SDK for PowerPoint Document
Convert. Convert PowerPoint to PDF. Extract, copy and paste PowerPoint Pages. Annotation & Thumbnail. Add and burn annotation to PowerPoint.
www.rasteredge.com
Chapter 31
Numerical methods
Several functions are available to aid in the construction of special-purpose estimators: one group
of functions are used to maximize user-supplied functions via numerical methods: BFGS, Newton–
Raphson andSimulatedAnnealing. Another relevant function is fdjac, which produces a forward-
difference approximation to the Jacobian.
31.1 BFGS
The BFGSmax function has two required arguments: a vector holding the initial values of a set of
parameters, and a call to a function that calculates the (scalar) criterion to be maximized, given
the current parameter values and any other relevant data. If the object is in fact minimization, this
function should return the negativeofthe criterion. On successful completion,BFGSmax returns the
maximized value of the criterion and the vector given via the first argument holds the parameter
values which produce themaximum. It is assumedhere that the objective function is a user-defined
function (see Chapter13) with the following general set-up:
function scalar ObjFunc (matrix *theta, matrix *X)
scalar val = ... # do some computation
return val
end function
The first argument contains the parameter vector and the second may be used to hold “extra”
values that are necessary to compute the objective function, but are not the variables of the opti-
mization problem. For example, if the objective function were a loglikelihood, the first argument
would contain the parameters and the second one the data. Or, for more economic-theory inclined
readers, if the objective function were the utility of a consumer, the first argument might contain
the quantities of goods and the second one their prices and disposable income.
Example 31.1: Finding theminimum of theRosenbrock function
function scalar Rosenbrock(matrix *param)
scalar x = param[1]
scalar y = param[2]
return -(1-x)^2 - 100 * (y - x^2)^2
end function
matrix theta = { 0, 0 }
set max_verbose 1
M = BFGSmax(&theta, Rosenbrock(&theta))
print theta
The operation of BFGS can be adjusted using the set variables bfgs_maxiter and bfgs_toler
275
Chapter 31. Numerical methods
276
(see Chapter21). In addition you can provoke verbose output from the maximizer by assigning a
positive value to max_verbose, again via the set command.
The Rosenbrock function is often used as a test problem for optimization algorithms. It is also
known as “Rosenbrock’s Valley” or “Rosenbrock’s Banana Function”, on account of the fact that its
contour lines are banana-shaped. It is defined by:
fx;y 1  x
2
100y  x
2
2
The function has a global minimum at x;y  1;1 where fx;y  0. Example31.1 shows a
gretl script that discovers the minimum usingBFGSmax (givinga verbose account ofprogress). Note
that, in this particular case, the function to be maximized only depends on the parameters, so the
second parameter is omitted from the definition of the function Rosenbrock.
Limited-memory variant
SeeByrdetal.(1995) (...FIXME: expand a little here ...)
Supplying analytical derivatives for BFGS
An optional third argument to the BFGSmax function enables the user to supply analytical deriva-
tives of the criterion function with respect to the parameters (without which a numerical approxi-
mation to the gradient is computed). This argument is similar to the second one in that it specifies
afunction call. In this case the function that is called must have the following signature.
Its first argument should be a pre-defined matrix correctly dimensioned to hold the gradient; that
is, if the parameter vector contains k elements, the gradient matrix must also be a k-vector. This
matrix argument must be given in “pointer” form so that its content can be modified by the func-
tion. (Note that unlike the parameter vector, where the choice of initial values can be important,
the initial values given to the gradient are immaterial and do not affect the results.)
In addition the gradient function must have as one ofits argument the parameter vector. This may
be given in pointer form (which enhances efficiency) but that isnot required. Additional arguments
may be specified if necessary.
Given the current parameter values, the function call must fill out thegradient vector appropriately.
It is not required that the gradient function returns any value directly; if it does, that value is
ignored.
Example31.2 illustrates, showing how the Rosenbrock script can be modified to use analytical
derivatives. (Note that since this is a minimization problem the values written into g[1] and g[2]
in the function Rosen_grad are in fact the derivatives of the negative of the Rosenbrock function.)
31.2 Newton–Raphson
BFGS, discussedabove, isan excellent all-purpose maximizer, andabout as robust as possible given
the limitations of digital computer arithmetic. The Newton–Raphson maximizer is not as robust,
but may converge much faster than BFGS for problems where the maximand is reasonably well
behaved—in particular, where it is anything like quadratic (see below). The case for usingNewton–
Raphson is enhanced if it is possible to supply a function to calculate the Hessian analytically.
The gretl function NRmax, which implements the Newton–Raphson method, hasa maximum offour
arguments. The first two (required) arguments are exactly as for BFGS: an initial parameter vector,
and a function call which returns the maximand given the parameters. The (optional) third argu-
ment is again asin BFGS:a function call that calculatesthe gradient. Specificto NRmax isan optional
fourth argument, namely a function call to calculate the (negative) Hessian. The first argument of
this function must be a pre-defined matrix of the right dimension to hold the Hessian—that is, a
kk matrix, where k is the length of the parameter vector—given in “pointer” form. The second
argument should be the parameter vector (optionally in pointer form). Other data may be passed
Chapter 31. Numerical methods
277
Example 31.2: Rosenbrock functionwithanalytical gradient
function scalar Rosenbrock (matrix *param)
scalar x = param[1]
scalar y = param[2]
return -(1-x)^2 - 100 * (y - x^2)^2
end function
function void Rosen_grad (matrix *g, matrix *param)
scalar x = param[1]
scalar y = param[2]
g[1] = 2*(1-x) + 2*x*(200*(y-x^2))
g[2] = -200*(y - x^2)
end function
matrix theta = { 0, 0 }
matrix grad = { 0, 0 }
set max_verbose 1
M = BFGSmax(&theta, Rosenbrock(&theta), Rosen_grad(&grad, &theta))
print theta
print grad
as additional arguments as needed. Similarly to the case with the gradient, if the fourth argument
to NRmax is omitted then a numerical approximation to the Hessian is constructed.
What is ultimately required in Newton–Raphson is the negative inverse of the Hessian. Note that
if you give the optional fourth argument, your function should compute the negative Hessian, but
should not invert it; NRmax takes care of inversion, with special handling for the case where the
matrix is not negative definite, which can happen far from the maximum.
Script31.3 extends the Rosenbrock example, using NRmax with a function Rosen_hess to compute
the Hessian. The functions Rosenbrock and Rosen_grad are just the same as in Example31.2 and
are omittedfor brevity.
The idea behind Newton–Raphson is to exploit a quadratic approximation to the maximand, under
the assumption that it is concave. If this is true, the method is very effective. However, if the
algorithm happens to evaluate the function at a point where the Hessian is not negative definite,
thingsmay go wrong. Script31.4 exemplifiesthis by usinga normal density, which is concave in the
interval  1;1 andconvexelsewhere. Ifthe algorithm isstartedfrom within the interval everything
goes well and NR is (slightly) more effective than BFGS. If, however, the Hessian is positive at the
starting point BFGS converges with only little more difficulty, while Newton–Raphson fails.
31.3 Simulated Annealing
Simulated annealing—as implemented by the gretl function simann—is not a full-blown maxi-
mization methodin its own right, but can be a useful auxiliary tool in problems where convergence
depends sensitively on the initial values ofthe parameters. The idea is that you supply initial values
and the simulated annealing mechanism tries to improve on them via controlledrandomization.
The simann function takes up to three arguments. The first two (required) are the same as for
BFGSmax and NRmax: an initial parameter vector and a function that computes the maximand. The
optional third argument is a positive integer giving the maximum number of iterations, n, which
Chapter 31. Numerical methods
278
Example 31.3: Rosenbrock function viaNewton–Raphson
function void Rosen_hess (matrix *H, matrix *param)
scalar x = param[1]
scalar y = param[2]
H[1,1] = 2 - 400*y + 1200*x^2
H[1,2] = -400*x
H[2,1] = -400*x
H[2,2] = 200
end function
matrix theta = { 0, 0 }
matrix grad = { 0, 0 }
matrix H = zeros(2, 2)
set max_verbose 1
M = NRmax(&theta, Rosenbrock(&theta), Rosen_grad(&grad, &theta),
Rosen_hess(&H, &theta))
print theta
print grad
Example 31.4: Maximization of a Gaussian density
function scalar ND(matrix x)
scalar z = x[1]
return exp(-0.5*z*z)
end function
set max_verbose 1
x = {0.75}
A = BFGSmax(&x, ND(x))
x = {0.75}
A = NRmax(&x, ND(x))
x = {1.5}
A = BFGSmax(&x, ND(x))
x = {1.5}
A = NRmax(&x, ND(x))
Chapter 31. Numerical methods
279
defaults to 1024.
Starting from the specified point in the parameter space, for each of n iterations we select at
random a new point within a certain radius of the previous one and determine the value of the
criterion at the new point. If the criterion is higher we jump to the new point; otherwise, we jump
with probability P (and remain at the previous point with probability 1   P). As the iterations
proceed, the system gradually “cools”—that is, the radius of the random perturbation is reduced,
as is the probability of making a jump when the criterion fails to increase.
In the course of this procedure n 1 points in the parameter space are evaluated: call them 
i
;i 
0;:::;n, where 
0
is the initial value given by the user. Let 
denote the “best” point among
1
;:::;
n
(highest criterion value). The value written into the parameter vector on completion is
then 
if 
is better than 
0
, otherwise 
n
. In other words, failing an actual improvement in
the criterion, simann randomizes the starting point, which may be helpful in tricky optimization
problems.
Example31.5 shows simann at work as a helper for BFGSmax in finding the maximum of a bimodal
function. Unaided, BFGSmax requires 60 function evaluations and 55 evaluations of the gradient,
while after simulated annealing the maximum is found with 7 function evaluations and 6 evalua-
tions of the gradient.1
Example 31.5: BFGSwith initialization via Simulated Annealing
function scalar bimodal (matrix x, matrix A)
scalar ret = exp(-qform((x-1)’, A))
ret += 2*exp(-qform((x+4)’, A))
return ret
end function
set seed 12334
set max_verbose on
scalar k = 2
matrix A = 0.1 * I(k)
matrix x0 = {3; -5}
x = x0
u = BFGSmax(&x, bimodal(x, A))
print x
x = x0
u = simann(&x, bimodal(x, A), 1000)
print x
u = BFGSmax(&x, bimodal(x, A))
print x
31.4 Computing a Jacobian
Gretl offers the possibility of differentiating numerically a user-defined function via the fdjac
function.
This function takes two arguments: an n1 matrix holding initial parameter values and a function
call that calculates and returns an m 1 matrix, given the current parameter values and any other
1
Your mileage may vary: these figuresare somewhatcompiler- and machine-dependent.
Documents you may be interested
Documents you may be interested