kalman.Rnw

\documentclass[10pt]{beamer}
%\usepackage[MeX]{polski}
\usepackage[utf8]{inputenc} 
%\usepackage[polish]{babel}    
\usepackage{parskip}
\usepackage{latexsym,gensymb,amsmath,amssymb,amsthm}
\usepackage{graphicx}
\usepackage{url}
\usepackage{caption}
\usepackage{wrapfig}
\usepackage{multicol}
\usepackage[QX]{fontenc}
\usepackage{lmodern}

\newcommand{\dxpt}{\ensuremath{\hat{x}_{t|t-1}}}
\newcommand{\dxtt}{\ensuremath{\hat{x}_{t-1|t-1}}}
\newcommand{\dxt}{\ensuremath{\hat{x}_{t|t}}}

\newcommand{\dPpt}{\ensuremath{{P}_{t|t-1}}}
\newcommand{\dPtt}{\ensuremath{{P}_{t-1|t-1}}}
\newcommand{\dPt}{\ensuremath{{P}_{t|t}}}



\usepackage[sorting=none,backend=biber,bibencoding=utf8]{biblatex} % load the package
\addbibresource{bibliografia.bib} % add a bib-reference file



\author{Mateusz Stachnik, Zygmunt Zawadzki}

\institute{Otwarte seminarium Katedry Statystyki}


\title{Dylematy związane ze stosowaniem Filtru Kalmana w prognozowaniu zjawisk ekonomicznych}
\begin{document}
\frame{
\titlepage
}


\section{Wprowadzenie}
\frame{\frametitle{Wprowadzenie} 
\tableofcontents
}


\frame{\frametitle{Materiały dodatkowe} 
Prezentacja powstała z wykorzystaniem R, Matlaba, \LaTeX 'a i knitr'a \footnote{\url{http://yihui.name/knitr/}}. Na uwagę zasługuje szczególnie knitr, który pozwolił w łatwy sposób zintegrować R i \LaTeX 'a.

Kod R użyty do wygenerowania przykładów w prezentacji można znaleźć w internecie pod adresem \url{https://github.com/zzawadz/KalmanSem}.

<<cache=FALSE,echo=FALSE, message=FALSE>>=
require(knitr)
opts_chunk$set(cache=TRUE,message=FALSE, warning=FALSE, echo = FALSE)
@

}

\section{Rys historyczny}


\frame{\frametitle{Rudolf Kalman} 

\begin{figure}
\includegraphics[scale = 0.4]{obrazki/Rudolf_Kalman}
\caption{Źródło: Wikipedia}
\end{figure}

\begin{itemize}
\item Urodzony w 1930 na Węgrzech. 
\item 1960 publikacja Filtru Kalmana.
\item Do dzisiaj wykłada.
\end{itemize}


}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Przestrzeń stanów}

\frame{\frametitle{Przestrzeń stanów - wprowadzenie - rys historyczny} 

W klasycznej teorii sterowania można było sterować układami o jednym wejściu i jednym wyjściu. Wraz z programem Apollo pojawiała się potrzeba sterowania wieloma obiektami jednocześnie. Problem polegał na tym, że występowały oddziaływania między obiektami.


%%% dodać np
\begin{align*}
Y = G \cdot X \\
G = \frac{1}{s+1}
\end{align*}
gdzie $G$ to transmitancja operatorowa układu, $Y$ to wyjście, $X$ wejście, $s$ - operator Laplace'a.
}

\frame{\frametitle{Przestrzeń stanów - wprowadzenie} 
Problem pojawiał się w opisie prostego układu o dwóch wejściach i dwóch wyjściach.
\begin{align*}
Y_1 = G_{11} \cdot X_1 + G_{12} \cdot X_2 \\
Y_2 = G_{21} \cdot X_1 + G_{22} \cdot X_2
\end{align*}

W podejściu przestrzeni stanów rozdziela się układ na dwie części. Część stanu, oraz część wyjścia. Powyższy układ można dzięki temu zapisać dwoma równaniami.
\begin{align*}
&\dot{x}(t) = A_cx(t) + B_cu(t) \\
&y(t) = H_cx(t) + D_cu(t)
\end{align*}
gdzie $A_c$ - macierz stanu, $B_c$ - macierz sterowania, $H_c$ - macierz wyjścia, $D_c$ - macierz transmisji. 

Oczywiste wydaje się, że zapis w przestrzeni stanów jest dużo prostszy i łatwo można go wykorzystać w dużych systemach. W dodatku jest bardziej ogólny niż zapis transmitancyjny. 

}


\frame{\frametitle{Przestrzeń stanów - układ dyskretny} 

Model dla czasu ciągłego postaci:
\begin{align*}
&\dot{x}(t) = A_cx(t) + B_cu(t) \\
&y(t) = H_cx(t) + D_cu(t)
\end{align*}
można również przedstawić w dziedzinie dyskretnej:
\begin{align*}
&x(t+1) = Ax(t) + Bu(t) \\
&y(t) = Hx(t)
\end{align*}

Istnieją metody pozwalające na przejście z czasu ciągłego do dyskretnego i z powrotem.

}

\frame{\frametitle{Przykład - estymacja pozycji pojazdu} 
 
Równania opisujące dynamikę układu:
\begin{align*}
v_{t+1} &= v_t + a_t T_p + w_v \\
s_{t+1} &= s_t + v_tT_p + \frac{1}{2}a_tT^2_p +w_s
\end{align*}

Układ zapisany w przestrzeni stanów:
\begin{align*}
\begin{bmatrix} s_{t+1} \\ v_{t+1} \end{bmatrix} &= \begin{bmatrix} 1 & T_p \\ 0 & 1 \end{bmatrix} \begin{bmatrix} s_{t} \\ v_{t} \end{bmatrix} + \begin{bmatrix} w_{st} \\ w_{vt} \end{bmatrix} \\
y_t &= \begin{bmatrix} 1 & 0\end{bmatrix} \begin{bmatrix} s_t \\ v_t \end{bmatrix} + z_t
\end{align*}
$T_p$ - czas próbkowania, równy $0.05s$

}


\frame{\frametitle{Przykład - estymacja pozycji pojazdu} 
\begin{center}
$\sigma_a = 0.1$,~~  $\sigma_s = 0.1$
\end{center}
\begin{figure}
\includegraphics[scale = 0.55]{obrazki/przegiegS}
\end{figure}
}


\frame{\frametitle{Cel filtracji} 

Głównym celem filtracji jest poznanie rzeczywistej wartości sygnału. Czasem problemem jest też fakt, że sygnał mocno zaszumiony jest trudny do przetwarzania, szczególnie tyczy się to różniczkowania takiego sygnału.
\begin{figure}
\includegraphics[scale = 0.6]{obrazki/filtr}
\end{figure}

}


\frame{\frametitle{Filtr analogowy} 

Filtry analogowe projektuje się zazwyczaj w domenie częstotliwości. 

\begin{figure}
\includegraphics[scale = 0.43]{obrazki/anal1}
\end{figure}

}

\frame{\frametitle{Filtr analogowy - dolnoprzepustowy drugiego rzędu} 
\begin{figure}
\includegraphics[scale = 0.43]{obrazki/anal2}
\end{figure}
}

\frame{\frametitle{Realizacja filtru dolnoprzepustowego jako filtru cyfrowego} 

\begin{align*}
G = \frac{1}{s^2 + s + 1}
\end{align*}
przechodząc z transmitancji ciągłej do dyskretnej metodą ZOH (ZERO ORDER HOLD) z czasem próbkowania 0.05s otrzymujemy:
\begin{align*}
G = \frac{0.0012z^{-1} + 0.0012 z^{-2}}{1 - 1.949z^{-1} + 0.095z^{-2}}
\end{align*}
Z transmitancji dyskretnej można w prosty sposób przejść na równanie różnicowe:
\begin{align*}
y_t = 0.0012x_{t-1} + 0.0012 x_{t-2} + 1.949y_{t-1} - 0.095y_{t-2}
\end{align*}

}


\frame{\frametitle{Charakterystyka czasowa - położenie} 
\begin{figure}
\includegraphics[scale = 0.6]{obrazki/anal3}
\end{figure}
}

\frame{\frametitle{Charakterystyka czasowa - prędkość} 
\begin{figure}
\includegraphics[scale = 0.6]{obrazki/anal4}
\end{figure}
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\frame{\frametitle{Filtr cyfrowy - średnia ruchoma} 

Prosta średnia ruchoma jest klasycznym filtrem cyfrowym typu MA:
\begin{align*}
y_t = \frac{1}{n}\sum^{t}_{i = t-n+1} x_{t-i}
\end{align*}

W rozważanym przypadku przyjęto $n = 10$.
 
}

\frame{\frametitle{Filtr cyfrowy - charakterystyka czasowa - położenie} 

\begin{figure}
\includegraphics[scale = 0.6]{obrazki/smooth1}
\end{figure}

}

\frame{\frametitle{Filtr cyfrowy - charakterystyka czasowa - prędkość} 
\begin{figure}
\includegraphics[scale = 0.6]{obrazki/smooth2}
\end{figure}
}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Filtr Kalmana}

\frame{\frametitle{Filtr Kalmana} 


Zadaniem filtracji Kalmana jest estymacja wartości wektora stanu x przy niepewnych parametrach obiektu na podstawie zmierzonych próbek obarczonych szumem.

Cały proces przebiega dwuetapowo - najpierw następuje faza predyckji na chwilę $t$, na podstawie chwili $t-1$, następnie na podstawie obserwacji (już w chwili t), następuje uaktualnienie estymowanych wartości.

Dla ustalenia uwagi:
\begin{equation}
\begin{split}
x_{t} = Ax_{t-1} + Bu_t + w_t \\
y_t = Hx_t + v_t
\end{split}
\end{equation}
\begin{tabular}{rr}

  $x_t$ - wektor zmiennych stanu   & $u_t$ - wektor sterowania - znany\\
  $y_t$ - wektor wyjścia - obserwowany & A,B,H - znane macierze \\
  $w_t \sim WN(0,W)$ & $v_t \sim WN(0,V)$ \\
  W i V - diagonalne macierze kowariancji. &

\end{tabular} 

}

\subsection{Faza predykcji}

\frame{\frametitle{Faza predykcji} 

Faza predykcji opisana jest następującymi wzorami:
\begin{equation}
\label{eq:predmean}
\dxpt = A\dxtt + Bu_t 
\end{equation}
\begin{equation}
\label{eq:predcov}
\dPpt = AP_{t-1|t-1}A' + W
\end{equation}
We wzorach pojawia się nowa macierz $\dPpt$ - jest to macierz kowariancji dla prognozy pod warunkiem chwili $t-1$, definiowana jako:
\begin{equation}
\label{eq:ppt1}
\begin{split}
\dPpt = E[(x_t-\dxpt)(x_t-\dxpt)']
\end{split}
\end{equation}
natomiast $\dPt$ to macierz kowariancji dla $x_t$ w chwili $t$ definiowana jako:
\begin{equation}
\label{eq:ppt1}
\begin{split}
\dPt = E[(x_t-\dxt)(x_t-\dxt)']
\end{split}
\end{equation}
W chwili uaktualniania $\dPtt$ jest znana - obliczana jest ona bowiem w chwili $t-1$.

}

\subsection{Faza uaktualnienia}

\frame{\frametitle{Faza uaktualnienia} 

W momencie dokonania pomiaru następuje uatktualnienie wartości zmiennych stanu, na podstawie następujących wzorów:

\begin{equation}
\label{eq:gain}
K_t = (\dPpt H')(H \dPpt H' + V)^{-1}
\end{equation}
\begin{equation}
\label{eq:mean}
\dxt = \dxpt + K_t(y_t - H\dxpt)
\end{equation}
\begin{equation}
\label{eq:cov}
\dPt = \dPpt - K_tH\dPpt 
\end{equation}
Najbardziej kluczową wartością jest $K_t$ nazywane wzmocnieniem Kalmana - określa wpływ zaobserwowanego odychenia wartości obserwowanej od predykcji ($y_t - H\dxpt$) na aktualizację wartości zmiennych stanu - im $K_t$ większe, tym większa waga będzie przykładana obserwacji.


}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\frame{\frametitle{Filtr Kalmana - charakterystyka czasowa - położenie} 

\begin{figure}
\includegraphics[scale = 0.6]{obrazki/kalman1}
\end{figure}

}

\frame{\frametitle{Filtr Kalmana - charakterystyka czasowa - prędkość} 
\begin{figure}
\includegraphics[scale = 0.6]{obrazki/kalman2}
\end{figure}
}

\section{Zestawienie filtrów}

\frame{\frametitle{Zestawienie filtrów- charakterystyka czasowa - położenie} 
\begin{figure}
\includegraphics[scale = 0.6]{obrazki/zestawienie1}
\end{figure}
}

\frame{\frametitle{Zestawienie filtrów- charakterystyka czasowa - prędkość} 
\begin{figure}
\includegraphics[scale = 0.6]{obrazki/zestawienie2}
\end{figure}
}

\frame{\frametitle{Zestawienie filtrów- błąd położenia} 
\begin{figure}
\includegraphics[scale = 0.6]{obrazki/zestawienie3}
\end{figure}
}

\frame{\frametitle{Zestawienie filtrów- błąd prędkości} 
\begin{figure}
\includegraphics[scale = 0.6]{obrazki/zestawienie4}
\end{figure}
}


\frame{\frametitle{Filtr Kalmana - współczynnik tłumienia szumów} 

\begin{align*}
WTS = \frac{RMS(x_R - x_P)}{RMS(x_R-x_F)} = \frac{RMS(e_P)}{RMS(e_F)}
\end{align*}
gdzie $x_R$ - pozycja rzeczywista, $x_P$ - pozycja zmierzona obarczona błędem, $x_F$ - pozycja po filtracji.

W związku z czym im większy współczynnik WTS tym lepsze tłumienie szumów pomiarowych i lepsza estymacja zmiennych stanu. 

}

\frame{\frametitle{Zestawienie filtrów - współczynnik tłumienia szumów - położenie} 
\begin{figure}
\includegraphics[scale = 0.6]{obrazki/WTS}
\end{figure}
}

\frame{\frametitle{Zestawienie filtrów - współczynnik tłumienia szumów - prędkość} 
\begin{figure}
\includegraphics[scale = 0.6]{obrazki/WTS2}
\end{figure}
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Przykład ekonomiczny}

\frame{\frametitle{Prosty przykład ekonomiczny - trend liniowy} 

Bardzo prostym i intuicyjnym przykładem jest model trendu liniowego - który możemy zapisać w przestrzeni stanów w postaci:
\begin{equation}
\begin{split}
\begin{bmatrix} \alpha \\ \beta \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} \alpha \\ \beta \end{bmatrix}\\
y(t) = \begin{bmatrix} 1 & 0 \end{bmatrix} \begin{bmatrix} \alpha \\ \beta \end{bmatrix} + v(t)
\end{split}
\end{equation}

Parametry $\alpha$ i $\beta$ są nieobserwowanymi zmiennymi stanu, które w procesie filtracji chcemy estymować.

}

\frame{\frametitle{Filtr - główna idea. Trend liniowy - równania.} 
\begin{align*}
&\begin{bmatrix} \alpha \\ \beta \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} \alpha \\ \beta \end{bmatrix}\\
&y(t) = \begin{bmatrix} 1 & 0 \end{bmatrix} \begin{bmatrix} \alpha \\ \beta \end{bmatrix} + v_t
\end{align*}
\begin{align*}
t = 1 \\
&y_1 = \alpha_1 + v_1 \\
&\alpha_1 = \alpha_0 + \beta_0 \\
&\beta_1 = \beta_0 \\
\\
t = 2 \\
&y_2 = \alpha_2 + v_2 = \alpha_0 + 2 \beta_0 \\
&\alpha_2 = \alpha_1 + \beta_1 = \alpha_0 + \beta_0 + \beta_1 = \alpha_0 + 2 \beta_0 \\
&\beta_2 = \beta_1 = \beta_0
\end{align*}

\textbf{Jaki więc zysk ma ekonomista z wykorzystania przestrzeni stanów?}

}


\frame{\frametitle{Model trendu liniowego ze zmieniającymi się parametrami} 

Co więcej przedstawienie modelu w przestrzeni stanów pozwala nam rozszerzyć klasę omawianych modeli. Zaprezentowany przepadek trendu możemy rozszerzyć na przypadek zmieniających się w czasie parametrów np. $\alpha$ i $\beta$ podlegają błądzeniu losowemu - wtedy równania stanu będą wyglądały następująco:

\begin{equation}
\begin{split}
\begin{bmatrix} \alpha_{t+1} \\ \beta_{t+1} \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} \alpha_t \\ \beta_t \end{bmatrix} + w_t\\
y_t = \begin{bmatrix} 1 & 0 \end{bmatrix} \begin{bmatrix} \alpha_t \\ \beta_t \end{bmatrix} + v_t
\end{split}
\end{equation}


}



\frame{\frametitle{Model trendu liniowego ze zmieniającymi się parametrami} 
<<echo=FALSE,fig.height = 3, fig.width = 4, message = FALSE>>=
source("KodR/exTrend.R")
@

}

\frame{\frametitle{Model trendu liniowego ze zmieniającymi się parametrami - realizacje i wartości estymowane} 
<<echo=FALSE,fig.height = 3, fig.width = 4, message = FALSE>>=
source("KodR/exTrendEstPar.R")
@

}


\frame{\frametitle{Model trendu liniowego - przykład empiryczny} 
<<echo=FALSE,fig.height = 3, fig.width = 4, message = FALSE>>=
source("KodR/estModelSP500.R")
@

}


\frame{\frametitle{Model trendu liniowego - przykład empiryczny - predykcja} 
<<echo=FALSE,fig.height = 4, fig.width = 6, message = FALSE>>=
source("KodR/estForecastSP500.R")
@

}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





\section{Dylematy?}

\frame{\frametitle{Dylematy?} 

\begin{itemize}
\item W modelach znanych z teorii sterowania większość parametrów przestrzeni stanu jest znana - są to własności wynikające z fizyki.
\item W przypadku modeli ekonomicznych w zasadzie wszystkie parametry podlegają estymacji, w tym rozmiar przestrzeni stanu (np. klasyczny problem związany z doborem zmiennych).
\end{itemize}

}

\section{Bibliografia}
\frame{\frametitle{Bibliografia - podstawy} 
\begin{refsection}

\nocite{commandeur2007introduction}
\nocite{kalie}
\nocite{konar}
\printbibliography[heading=subbibliography]

\end{refsection}
}


\frame{\frametitle{Bibliografia - układy dynamiczne} 
\begin{refsection}

\nocite{ogata}
\nocite{osowski}
\printbibliography[heading=subbibliography]

\end{refsection}
}



\frame{\frametitle{Bibliografia - Filtr Kalmana w R} 
\begin{refsection}
\nocite{Tusell}
\nocite{Petris2}
\nocite{Petris}
\printbibliography[heading=subbibliography]
\end{refsection}
}




\end{document}