Library: | kalman |
See also: | ksmoother kem kfilter kfilter2 calibrIC ICerzsep ICerz rlsfil kemitor kemitor2 |
Quantlet: | rICfil | |
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} = rICfil(y,mu,Sig,H,F,Q,R,cliptyp,A,b) | |
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 | |
cliptyp | numeric : either 0 for simultaneous clipping 1 for separate clipping (AO) 2 for separate clipping (IO) | |
As | T,n,n matrix, "vector" of Lagrange-Matrices A; are to be produced by calibrIC or n,n matrix, uniform Lagrange-Matrix A | |
bs | T vector, the clipping heights; are to be produced by calibrIC 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 fact=1.2 expl=2 A0=0 ;; produces default values inside the called quantlets b0=-1 ;; produces default values inside the called quantlets typ= 0 ;; simultaneous clipping ergIC=calibrIC(T,Sig,H,F,Q,R,typ,A0,b0,e,N,eps,itmax,expl,fact,aus) ;; gets the clipping heights b and the Lagrange Multiplyers A to efficiency loss e A=ergIC.A b=ergIC.b res = kfilter2(y,mu,Sig, H,F,Q,R) fx = res.filtX res= rICfil(y,mu,Sig,H,F,Q,R,typ,A,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, rIC-sim, 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=#(24,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=300 eps=0.01 itmax=15 aus=4;; lots of output;; may be set to 0 fact=1.2 expl=2 A0=0 ;; produces default values inside the called quantlets b0=-1 ;; produces default values inside the called quantlets typ= 0 ;; simultaneous clipping ergIC=calibrIC(T,Sig,H,F,Qid,R,typ,A0,b0,e,N,eps,itmax,expl,fact,aus) ;; gets the clipping heights b and the Lagrange Multiplyers A to efficiency loss e A=ergIC.A b=ergIC.b res = kfilter2(y,mu,Sig, H,F,Qid,R) fx = res.filtX res= rICfil(y,mu,Sig,H,F,Qid,R,typ,A,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, rIC-sim, 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=300 eps=0.01 itmax=15 aus=4;; lots of output;; may be set to 0 fact=1.2 expl=2 A0=0 ;; produces default values inside the called quantlets b0=-1 ;; produces default values inside the called quantlets typ= 1 ;; separate clipping(AO) ergIC=calibrIC(T,Sig,H,F,Q,R,typ,A0,b0,e,N,eps,itmax,expl,fact,aus) ;; gets the clipping heights b and the Lagrange Multiplyers A to efficiency loss e A=ergIC.A b=ergIC.b res= kfilter2(y,mu,Sig, H,F,Q,R) fx = res.filtX fy=(H*fx')' res= rICfil(y,mu,Sig,H,F,Q,R,typ,A,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.)