| Library: | kalman |
| See also: | ksmoother kem kfilter kfilter2 calibrLS rICfil kemitor kemitor2 |
| Quantlet: | rlsfil | |
| Description: | Calculates a filtered time serie (uni- or multivariate) using a robust, recursive Filter based on LS-optimality, the rLS-filter. additionally to the Kalman-Filter one needs to specify the degree of robustness one wants to achieve; this is done either by specifying a clipping height or by specifying a relative loss w.r.t. the classical Kalman Filter in the ideal model in terms of MSE. The state-space model is assumed to be in the following form: y_t = H x_t + v_t x_t = F x_t-1 + w_t x_0 ~ (mu,Sig), v_t ~ (0,Q), w_t ~ (0,R) All parameters are assumed known. |
| Usage: | {filtX,KG,PreviousPs,clipInd} = rlsfil(y,mu,Sig,H,F,Q,R,bs) | |
| Input: | ||
| y | T x m matrix of observed time series, T is the number of observations, m is the dimension of time series | |
| mu | n x 1 vector, the mean of the initial state | |
| Sig | n x n covariance matrix of the initial state | |
| H | m x n matrix | |
| F | n x n matrix | |
| Q | m x m variance-covariance matrix | |
| R | n x n variance-covariance matrix | |
| bs | T vector, the clipping heights ; are to be produced by calibrLS or numerical (uniform) clipping height | |
| Output: | ||
| filtX | T x n matrix of filtered time series | |
| KG | T x n x m "vector" of Kalman-Gain-matrices for different times | |
| PreviousPs | T x n x n "vector" of P_t|t's | |
| clippInd | T vector of 1 (if rlsfil clipped at instance t) and 0 (else) | |
library("kalman")
library("plot")
serie = read("kalmanao.dat")
; produced from kalman1.dat by manipulating observations 50, 60, 90
y = serie[,2]
mu = 10
Sig = 0
H = 1
F = 1
Q = 9
R = 9
T=dim(y)
e=0.05
N=100
eps=0.01
itmax=15
aus=4 ;lots of output, may be set to 0
ergLS=calibrLS(T,Sig,H,F,Q,R,e,N,eps,itmax,aus)
; gets the clipping heights to efficiency loss e
b=ergLS.b
res = kfilter2(y,mu,Sig, H,F,Q,R)
fx = res.filtX
res= rlsfil(y,mu,Sig, H,F,Q,R,b)
frx = res.filtX
origy= serie[,1]~serie[,2]
origy= setmask(origy, "line", "blue", "thin")
fx = serie[,1]~fx
fx = setmask(fx, "line", "red", "thin")
frx = serie[,1]~frx
frx = setmask(frx, "line", "green", "thin")
clip=serie[,1]~(res.clipInd)
clip=paf(clip,clip[,2]==1)
clip[,2]=0
setmaskp(clip,4, 3, 4)
disp = createdisplay(1,1)
show(disp,1,1, origy,fx,frx,clip)
setgopt(disp,1,1, "title", "KalmanData1 + AO's in t=50,60,90")
setgopt(disp,1,1, "xlabel", "t")
setgopt(disp,1,1, "ylabel", "y, rLS, Kalman")
Serie of observations is displayed with blue colour, Kalman-filtered serie is displayed with red colour, rls-filtered serie is displayed with green colour. The stages rlsfil has clipped at are marked by a flag.
library("kalman")
library("plot")
T=50
e=0.1
mu = #(20,0)
Sig = #(0,0)~#(0,0)
H = #(0.3,-0.3)~#(1,1)
F = #(1,0)~#(1,0)
R = #(0,0)~#(0,9)
mid=#(0,0)
Qid= #(9,0)~#(0,9)
mcont=#(25,30)
Qcont=0.1.*Qid
AOr=0.1
randomize(0)
ErrX = epscontnorm(T,0,mid,R,mcont,Qcont,0) ;; no contamination in ErrX
ErrY = epscontnorm(T,AOr,mid,Qid,mcont,Qcont,0) ;; 10% contamination in ErrY
sim = kemitor2(T,mu,H,F,(ErrY.X),(ErrX.X));; simulates a length 100 path of this state space model under AO's
y=sim.Y ;; observation
Xact=sim.X ;; actual state
N=10000
eps=0.01
itmax=15
aus=4 ;lots of output, may be set to 0
ergLS=calibrLS(T,Sig,H,F,Qid,R,e,N,eps,itmax,aus)
; gets the clipping heights to efficiency loss e
b=ergLS.b
res = kfilter2(y,mu,Sig, H,F,Qid,R)
fx = res.filtX
res= rlsfil(y,mu,Sig, H,F,Qid,R,b)
frx = res.filtX
i=(1:T)
Xact1 = i~(Xact[,1])
Xact1 = setmask(Xact1, "line", "blue", "thin")
fx1 = i~(fx[,1])
fx1 = setmask(fx1, "line", "red", "thin")
frx1= i~(frx[,1])
frx1 = setmask(frx1, "line", "green", "thin")
ym1=max(vec((Xact[,1])~(fx[,1])~(frx[,1]) ) ) ;Ind-flags on top of graphics
ym2=min(vec((Xact[,1])~(fx[,1])~(frx[,1]) ) ) ;clip-flags on bottom of graphics
flagInd=i~(ErrY.Ind)
flagInd=paf(flagInd,flagInd[,2]==1)
flagInd[,2]=ym1*((ym1>0)*1.1+(ym1<0)*0.9)
flagclip=i~(res.clipInd)
flagclip=paf(flagclip,flagclip[,2]==1)
flagclip[,2]= ym2*((ym2<0)*1.1+(ym2>0)*0.9)
setmaskp(flagInd,4, 3, 4)
setmaskp(flagclip,2, 4, 4)
disp = createdisplay(1,2)
show(disp,1,1,Xact1,fx1,frx1,flagInd,flagclip)
setgopt(disp,1,1, "title", "simulated Model under AO -- 1st coord.")
setgopt(disp,1,1, "xlabel", "t")
setgopt(disp,1,1, "ylabel", "x, rLS, Kalman")
Xact2 = i~(Xact[,2])
Xact2 = setmask(Xact2, "line", "blue", "thin")
fx2 = i~(fx[,2])
fx2 = setmask(fx2, "line", "red", "thin")
frx2= i~(frx[,2])
frx2 = setmask(frx2, "line", "green", "thin")
ym1=max(vec((Xact[,2])~(fx[,2])~(frx[,2]) ) )
ym2=min(vec((Xact[,2])~(fx[,2])~(frx[,2]) ) )
flagInd[,2]=ym1*((ym1>0)*1.1+(ym1<0)*0.9)
flagclip[,2]= ym2*((ym2<0)*1.1+(ym2>0)*0.9)
setmaskp(flagInd,4, 3, 4)
setmaskp(flagclip,2, 4, 4)
show(disp,1,2,Xact2,fx2,frx2,flagInd,flagclip)
setgopt(disp,1,2, "title", "simulated Model under AO -- 2nd coord.")
setgopt(disp,1,2, "xlabel", "t")
setgopt(disp,1,2, "ylabel", "x, rLS, Kalman") ;
Original serie of states is displayed with blue colour, conventionally filtered serie is displayed with red colour. robustly filtered serie is displayed with green colour. Red Flags are displayed where actual AO's where simulated and green ones where the robustly filtered serie was clipped. (y is a lagged bivariate MA process with errors.)
library("kalman")
library("plot")
serie = read("kalman2.dat")
y = serie[,2]
mu = #(0,0)
Sig = #(0,0)~#(0,0)
H = #(1,0)'
F = #(0.5,1)~#(-0.3,0)
R = #(1,0)~#(0,0)
Q = 4
T = dim(y)
e=0.05
N=10000
eps=0.01
itmax=15
aus=4;; lots of output;; may be set to 0
ergLS=calibrLS(T,Sig,H,F,Q,R,e,N,eps,itmax,aus)
; gets the clipping heights to efficiency loss e
b=ergLS.b
res= kfilter2(y,mu,Sig, H,F,Q,R)
fx = res.filtX
fy=(H*fx')'
res= rlsfil(y,mu,Sig, H,F,Q,R,b)
frx = res.filtX
fry=(H*frx')'
origy = serie[,1]~serie[,2]
origy = setmask(origy, "line", "blue", "thin")
fy = serie[,1]~fy
fy = setmask(fy, "line", "red", "thin")
fry = serie[,1]~fry
fry = setmask(fry, "line", "green", "thin")
flags = serie[,1]~(res.clipInd)
flags=paf(flags,flags[,2]==1)
flags[,2]=0
setmaskp(flags,4, 3, 4)
disp = createdisplay(1,1)
show(disp,1,1,origy,fy,fry,flags)
setgopt(disp,1,1, "title", "KalmanData2 in Observation Space")
setgopt(disp,1,1, "xlabel", "t")
setgopt(disp,1,1, "ylabel", "y, y-rLS, y-Kalman")
Original serie of observations displayed with blue colour, kalman-filtered serie y_t|t is displayed with red colour, rLS-filtered serie y_t|t is displayed with green colour. Red Flags are displayed where actual AO's where simulated and where the robustly filtered serie was clipped. (y is an AR(2) process with errors.)