Code for Figure 5.16

main.m


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58%
% Show 95 percent confidence ellipses versus time
% jbr
% 3/30/2013

ntsteps = 10;
U=0.5*[sqrt(2), sqrt(2); sqrt(2), -sqrt(2)];

theta = (-25)*360/(2*pi);
A = [cos(theta), sin(theta); -sin(theta), cos(theta)];
n = rows(A);
C = [0.5 0.25];
G = eye(2);
m = 2;
B = eye(m);
g = columns(G);
p = rows(C);
alpha = 0.95;
nell = 281;
bn = chi2inv(alpha, n);

prior = 3;
% Q0 = prior*U*diag([3, 1/2])*U';
Q0 = prior*eye(n);
Q  = 0.01*eye(g);
R  = 0.01*eye(p);

%[Ls, Pmins, Ps, Es] = dlqe(A, G, C, Q, R);

% transient solution
Pminus(1) = {Q0};

for i = 1:ntsteps
L{i} = Pminus{i}*C'*inv(R+C*Pminus{i}*C');
P{i} = Pminus{i} - L{i}*C*Pminus{i};
[exm{i}(:,1), exm{i}(:,2)] = ellipse(inv(Pminus{i}),bn, nell);
[ex{i}(:,1), ex{i}(:,2)] = ellipse(inv(P{i}),bn, nell);
if (i == ntsteps) break end
Pminus{i+1} = A*P{i}*A' + G*Q*G';
end

figure()
hold on
for i = 1:ntsteps
plot (ex{i}(:,1),ex{i}(:,2))
end
axis ([-1,1,-2,2]);

figure()
hold on
for i = 1:ntsteps
plot (exm{i}(:,1),exm{i}(:,2))
end
axis();

exx = cell2mat(ex);
save KFpayoff.dat  exx