2015年1月15日 星期四

[線性系統] 離散時間系統的可觀察性質 (Observability)

考慮 離散時間 線性非時變 (Discrete Time Linear Time Invariant, DT-LTI)系統
\[\begin{array}{l}
x(k + 1) = Ax(k)\\
y(k) = Cx(k)
\end{array}\] $x \in \mathbb{R}^n, y \in \mathbb{R}^p$。

我們說上述系統為 可觀察 (observable) 或稱 $(A,C)$ 可觀察 若下列條件成立:
存在 常數 $N < \infty$ ,使得對任意 初始狀態 $x(0)$ 而言,可用 $N$ 組量測輸出 $\{y(0), y(1),...,y(N-1)\}$  uniquely 決定該初始狀態 $x(0)$。

Comment
1. 上述定義可類比 可控制性條件,

2. 事實上若我們無法透過 $n$ 組 量測輸出 來區別 $x(0)$ 則就算給額外再多的量測輸出 e.g., $N>n$ 組 仍無法區別 $x(0)$。(此結果可由 Cayley-Hamilton Theorem 證明。)

以下我們看個 unobservable 的例子
上方的方塊圖 顯示了 子系統 $G(z)$ 的狀態 無法從輸出 $Y(z)$ 觀察到。 (圖中的 $(z)$ 表示對原系統做 Z-transform。)



觀察性基本問題:
透過 sensor 所量測到的輸出 $y$ 是否足夠讓我們找出系統 初始狀態 $x(0)$ uniquely?

為何我們關心 初始狀態? 因為一但有初始狀態則其餘任意時刻狀態均可透過狀態方程求解獲得。亦即 給定 $x(0)$ 則
\[\left\{ \begin{array}{l}
x(1) = Ax(0)\\
x(2) = {A^2}x(0)\\
...\\
x\left( N \right) = {A^N}x\left( 0 \right)
\end{array} \right.\]故若給定初始狀態 $x(0)$ 則其餘任意時刻狀態 $x(1), x(2),...,x(N)$均可透過狀態方程 $x(k+1) = Ax(k)$ 獲得。

但現在我們僅給定 $y(0),...,y(N)$ 亦即我們僅知道
\[ \Rightarrow \left\{ \begin{array}{l}
y(0) = Cx(0)\\
y(1) = Cx(1) = CAx(0)\\
y(2) = Cx(2) = C{A^2}x(0)\\
...\\
y\left( N \right) = C{A^N}x\left( 0 \right)
\end{array} \right.\]或者更進一步改寫成矩陣形式
\[\left[ {\begin{array}{*{20}{c}}
{y(0)}\\
{y(1)}\\
 \vdots \\
{y(N)}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
C\\
{CA}\\
 \vdots \\
{C{A^N}}
\end{array}} \right]x\left( 0 \right)\] 想問是否可從這些 $\{y(0),...,y(N)\}$ 回推 $x(0)$ (uniquely!)。 如果可以我們稱系統可觀察,若不行我們稱系統不可觀察。

Comment:
回憶在線性代數中,我們說 $Ax = b$ 解存在 若且唯若 $A$ 有 indepenent row (此對應 controllabilility problem);若我們說 $Ax = b$ 有唯一解 (注意 唯一不保證存在!!),若且為若 $A$ 有 independent column (此對應 observability problem)。

故我們要求觀察性矩陣 $O$ ( 其維度 $\dim(O) = Np  \times n$)
\[O = \left[ {\begin{array}{*{20}{c}}
C\\
{CA}\\
 \vdots \\
{C{A^N}}
\end{array}} \right]\] 有 independent column。由 Caley-Hamilton Theorem ,我們可僅考慮 $n$ 個量測輸出,則我們有
\[\left[ {\begin{array}{*{20}{c}}
{y(0)}\\
{y(1)}\\
 \vdots \\
{y(n - 1)}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
C\\
{CA}\\
 \vdots \\
{C{A^{n - 1}}}
\end{array}} \right]x\left( 0 \right)\]故觀察性矩陣 $O$ ( 其維度 $\dim(O) = np  \times n$)
\[O = \left[ {\begin{array}{*{20}{c}}
C\\
{CA}\\
 \vdots \\
{C{A^{n - 1}}}
\end{array}} \right]\]需要有 full column rank $n$ 我們總結 可觀察性的測試如下:

Theorem: Observability Rank Test
若系統 $(A,C)$ 為 observable 若且唯若 $rank(O) = n$。

同樣的我們也有 Hautus Lemma for observability

Lemma: Hautus Lemma for observability
一個系統為 observable 若且唯若 對任意 $\lambda \in \mathbb{C}$,
\[rank\left[ {\begin{array}{*{20}{c}}
{\lambda I - A}\\
C
\end{array}} \right] = n\]
注意到上述 Lemma 中,若 $\lambda \notin eig(A)$,則前面 $n$ rows 為 linearly independent ,故我們可已不用檢驗整個複數平面 $\lambda \in \mathbb{C}$,僅僅需要檢驗 $\lambda \in eig(A)$ 的部分即可。故我們得到以下修正引理:

Lemma: Modified Hautus Lemma
一個系統為 observable 若且唯若 對任意 $\lambda \in eig(A)$,
\[rank\left[ {\begin{array}{*{20}{c}}
{\lambda I - A}\\
C
\end{array}} \right] = n\]







2015年1月8日 星期四

[博弈論] Kelly criterion - Simplest Case

這次要介紹 凱利 J. Kelly, Jr 在 1956 年利用 Information Theory 所提出的一套 賭博理論結果,又稱作 凱利判準 (Kelly Criterion) ,此結果其後又被其在同僚 E. D. Throp 進一步推廣,且在近幾十年中已被逐步推廣且應用在 避險基金 與 投信銀行 等金融業中。以下我們將簡介最簡單的 Kelly criterion 型式:

凱利判準 (Kelly Criterion): 賭徒參與賭局,應追求 最大化 長期 報酬增長率 $G$ (asymptotically maximize the growth rate of wealth)
凱利公式 (Kelly formula): 此公式用來計算 每次賭金押注應該是多少可達成最大化 長期報酬增長率 $G$。

不過在介紹之前我們先考慮以下一個簡單的賭局情形:

Example
考慮 某賭徒身懷 $V_0$ 元全身資產 興致勃勃的 參與 某 投擲銅板 賭局,若銅板出現正面,則賭徒可贏得 $1$ 元,反之若出現反面則 賭徒輸掉 $1$元。且每次投擲銅板之間彼此互為獨立,現在給定 銅板出現正面機率為 $p$ 則反面出現的機率為 $1-p$ 並且 定義隨機變數 $X_k$表示 第 $k$ 次 投擲銅板,若銅板出現正面 我們記做 $X_k = 1$ 反之則記做 $X_k = -1$。

試問 (a) 賭徒在 第 $k$ 次 投擲銅板 的獲利期望值為何 $E[X_k]=?$
(b) 令第 $k$ 次賭注為 $B_k$ $(k=1,...,n)$,試問累計到 $n$ 次投擲銅板時 的獲利期望值為何 $E[V_n]=?$

Solution
(a) 第 $k$ 次 投擲銅板 的獲利期望值 由定義可知
\[\begin{array}{l}
E\left[ {{X_k}} \right] = 1 \cdot P\left( {{X_k} = 1} \right) + \left( { - 1} \right) \cdot P\left( {{X_k} =  - 1} \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}\;
\end{array} = 1 \cdot p + \left( { - 1} \right) \cdot \left( {1 - p} \right) = 2p - 1
\end{array}\]
(b) 對 $k=1,2,...$ 我們有 $V_k = V_{k-1} + X_kB_k $ 故
\[{V_n} = {V_0} + \sum\limits_{k = 1}^n {B_k}{X_k}
\]現在對上式取期望值可得
\[\begin{array}{l}
E{V_n} = {V_0} + \sum\limits_{k = 1}^n E [{B_k}{X_k}]\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = {V_0} + \sum\limits_{k = 1}^n {E[{B_k}]E[{X_k}]} \ \ \ \ \ (*) \\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = {V_0} + \sum\limits_{k = 1}^n {\left( {2p - 1} \right)E[{B_k}]}
\end{array}\]

Comments:
1. 注意到若 $2 p -1 >0$ 則表示第 $k$ 次 投擲銅板賭局的 獲利的期望值為正 ($EX_k >0$)。
2. 式 $(*)$ 中利用了 $B_k$ 與 $X_k$ 的獨立性 (儘管 $B_k$ 與 $X_{k-1}$ 有關 )。


有了以上想法我們開始檢驗以下幾種策略。假設你已知此投擲銅板的賭局的一個 "內線"消息,也就是此銅板是不公平的銅板,其出現正面的機率為 $p > \frac{1}{2}$ ,若你是賭徒試問應如何建構必勝法?

===============
策略 A: 既然已知出現正面機率較高,應該一次就押注全部的錢就對了。
===============
ANS: 若採用此策略,不難想像如果這位賭徒運氣不佳,前面好幾次都連續出現 反面,則 賭徒如果採用此全部押注的策略,只要出現一次反面就會輸個精光。且如果你說運氣很好 比如說此賭徒連贏 $n$ 次機率 為 $p^n $ 且 $0 < p < 1$ 故現在若 $n \to \infty$ 可想而知 賭徒連續贏無窮次的機率為 $0$ almost surely,亦即此賭徒必定破產。

對於 策略A 的結論如下:
只要輸一次賭徒就破產,且長期而言 ($n \to \infty$ ) 賭徒連贏的機率是 $0$,故若將全部資金投入賭局絕非必勝法。


故此很明顯我們該採取 將資金分批賭,那麼該如何分配這些資金才能 最大化賭徒贏錢的機會? 或者說 創造某種必勝法呢?

===============
策略 B: Martingale ! 每輸一次就加倍賭金。
===============
ANS: 此法看似必勝法事實上要求賭徒具備無窮資本。相關細節請讀者參閱本 BLOG 相關 Martingale 理論的介紹。


===============
策略 C: 利用 Kelly criterion 幫助我們決定每次賭注金的金額,使得長期而言,賭徒可贏得最大化 Geometric mean 收益。
===============

Kelly 定義 賭徒 參與第 $N$ 次賭局的資金 長期成長率 $G$ 為
\[G: = \mathop {\lim }\limits_{N \to \infty } \frac{1}{N}{\log _2}\left( {\frac{{{V_N}}}{{{V_0}}}} \right)\] 其中 $V_N$ 為 $N$ 次賭局後賭徒的手上的資金,$V_0$ 為 賭徒的初始資金。(NOTE: 我們在此將成長率定為 $\log_2$ 是依照 Kelly 1956 原文,事實上亦可直接訂為 $\ln(\cdot)$)

由於賭徒事先知道一個 "內線"消息,也就是此銅板是不公平的銅板,其出現正面的機率為 $p > \frac{1}{2}$ ,且此時賭徒只考慮投注比率為 $K$ ,而 $W$ 與 $L$ 各自代表贏得賭局 或者 輸掉賭局的次數,那麼考慮  $N$ 次賭局之後 (注意 $W + L = N$),賭徒剩餘資金可表示為
\[
V_N = (1 + K)^W(1 - K)^L V_0
\]且賭徒資金成長率 $G$ 亦可計算如下
\[\begin{array}{l}
G = \mathop {\lim }\limits_{N \to \infty } \frac{1}{N}{\log _2}\left( {\frac{{{V_N}}}{{{V_0}}}} \right)\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = \mathop {\lim }\limits_{N \to \infty } \frac{1}{N}{\log _2}\left( {\frac{{{{(1 + K)}^W}{{(1 - K)}^L}{V_0}}}{{{V_0}}}} \right)\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = \mathop {\lim }\limits_{N \to \infty } \frac{1}{N}{\log _2}\left( {{{(1 + K)}^W}{{(1 - K)}^L}} \right)\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = \mathop {\lim }\limits_{N \to \infty } \frac{1}{N}\left[ {W{{\log }_2}(1 + K) + L{{\log }_2}(1 - K)} \right]\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = \mathop {\lim }\limits_{N \to \infty } \left[ {\frac{W}{N}{{\log }_2}(1 + K) + \frac{L}{N}{{\log }_2}(1 - K)} \right]\\
\begin{array}{*{20}{c}}
{}&{}
\end{array} = p{\log _2}(1 + K) + \left( {1 - p} \right){\log _2}(1 - K)
\end{array}\]注意到上式 $-G$ 為 convex function of $K$,故可對其求解最佳 $K$ 值 (透過一階必要條件 $dG/dK =0$ )可得
\[\begin{array}{l}
\frac{{dG}}{{dK}} = 0\\
 \Rightarrow  - \frac{{1 - p}}{{\left( {1 - K} \right)\log \left[ 2 \right]}} + \frac{p}{{\left( {1 + K} \right)\log \left[ 2 \right]}} = 0\\
 \Rightarrow K = 2p - 1
\end{array}\]亦即 Kelly criterion 指出 若賭徒每次都以 $K = 2p-1$ 比率的賭金進行賭局,則可以最大化 Geometric mean 報酬率。且賭金成長率的最大值 $G_{max}$ 如下
\[\begin{array}{l}
{G_{\max }} = {\left. G \right|_{K = 2p - 1}} = p{\log _2}(2p) + \left( {1 - p} \right){\log _2}(2 - 2p)\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = p\left[ {{{\log }_2}(2) + {{\log }_2}(p)} \right] + \left( {1 - p} \right)\left[ {{{\log }_2}2 + {{\log }_2}(1 - p)} \right]\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}
\end{array} = 1 + p{\log _2}(p) + \left( {1 - p} \right){\log _2}(1 - p)
\end{array}\]

故以前例而言,若已知 $p = 0.55$ (因為已知小道消息幫助我們確認正面出現機率較高) 且賭徒身上帶著 $V_0$ 元作為賭本, 則 每次賭注金 應為 $(2 p - 1) V_0 = 0.1 V_0$ 亦即每次下注 $10 \%$。
Comments
以前述討論為例,若賭徒每次下注都低於 $10 \%$,則長期而言賭徒仍會持續勝利,但是獲利的成長率會較 每次都下注 $10 \%$ 來的慢 (因為非最佳解)。


上述的結果可以進一步推廣如下:首先引入 賠率 (賠率 := 贏的金額 / 輸的金額 );比如說下注金 $2$元,若賭贏則贏得 $4$ 元,賭輸則輸掉 $2$ 元。則此時賠率為 $4/2 = 2$ 。

Generalization of Kelly formula for Uneven payoff game  (1984, Thorp)
若已知賠率 $B$ 與獲勝機率 $p >0$ 且 $(B+1)p-1 >0$,則前述的 Kelly formula 可修正為
\[
K = \frac{(B + 1)p - 1} {B}
\] 上述結果來自於對下式最大化:
\[
\max_K G(K)
\]其中 $G(K) = p\ln \left( {1 + BK} \right) + \left( {1 - p} \right)\ln \left( {1 - K} \right)$

Example: 考慮賠率為 $B = 2$ 且 $p = 0.55$ 則 ,每次的最佳下注金比率應為
\[K = \frac{{(B + 1)p - 1}}{B} = \frac{{(2 + 1)\left( {0.55} \right) - 1}}{2} = 0.325
\]

Comment:
1. 事實上,讀者可以在文獻中找到更保守的投資方法,稱為 Half-Kelly (or fractiaonal-Kelly),一類常見的做法是限定 $ K^\star := 1/2 K^*$ 其中 $K^*$ 為 Kelly 最佳解。但不論是 Kelly 或者 Half-Kelly,此法則都還是有相當大的限制與問題:首先是 當賭徒把最佳解調降來得到 frctional-Kelly optimum 本身便相當具有爭議,此作法大概可以說成是把最佳化方法丟掉然後隨便亂調一個夠小的 Kelly optimum 並聲稱此解仍為最佳,既然如此何必費工做最佳化呢? 再者, Kelly Criterion 所求得的最佳比率 $K^*$ 儘管 "長期" 而言能最大化成長率,但不保證有限期間的績效,也就是說很有可能在有限期間透過 $K^*$ 去投資的 sampled path 仍會發生在某時刻大幅度資產減少 (亦即 會有極為巨大的 最大跌幅 (Max Drawdown) ),在 [4] 中我們證明就算是最基本的投擲銅板的情況,其 期望的最大跌幅仍超過 50% 以上。此現象在許多文獻中都被提及,一般將此性質稱為 Kelly optimum 具有過度投資的特性 (或稱 Kelly Criterion 是 too aggressive)。

2. 許多文獻試圖應用 Kelly Criterion 到股票市場之中,但實際應用上仍有非常多的限制,比如說有部分文獻透過 Talyor Approximation method 來求解最佳比率 $K^*$ 而無視 凱莉問題本質上是 concave program, 此類近似法可以給出漂亮的解析解並且看似合理,但其實仍潛藏諸多限制與謬誤。有興趣的讀者可參考個人的著作 [4].

3. 應用 Kelly Criterion 到股票市場有另外一類更嚴重的問題:回憶若 Kelly Criterion 應用在賭局情況,則賭徒大可假設 報酬的機率分佈已知,但是股票市場的機率分佈是完全未知,連最基本的 i.i.d. 或者 stationarity 報酬假設都僅只有在交易時間不長的情況下成立。有興趣的讀者在查閱相關文獻不難發現若貿然應用 Kelly Criterion 並透過模擬方式宣稱回測結果幾乎都無可避免上述的嚴重謬誤。


ref:
[1] Kelly, J. L., A New Interpretation of Information Rate," Bell System Technical Journal, 1956
[2] Thorp, E. O., "The Mathematics of Gambling," Lyle Stuart, Secaucus,NJ. 1984
[3] Thorp, E. O., The Kelly Criterion in Blackjack Sports Betting, and The Stock Market, 2007
[4] Hsieh, C.H., and Barmish, B. R., "On Kelly Betting: Some Limitations,"  Proceedings of 53rd Annual Allerton Conference, pp. 165-172, Monticello, IL., 2015

2015年1月4日 星期日

[隨機分析] Euler-Maruyama 法求 隨機微分方程 數值解 (利用 MATLAB)

在此我們介紹 Euler-Maruyama Method 求解 隨機微分方程   (Stochastic Differential Equation, SDE),現在 考慮  SDE 的積分型式 可寫為:
\[
X(t) = X_0 + \int_0^t f(X(s))ds + \int_0^t g(X(s))dW(s), \;\; 0 \le t \le T
\]其中 $f, g$ 為 純量函數 且 初始值 $X_0$ 為隨機變數;另外 第二項積分為 Ito integral。上式若有解則其解 $X(t)$ 對任意 $t$ 皆為隨機變數。

一般而言,上式可改寫為較為簡潔的 隨機微分方程的微分形式:
\[
dX(t) = f(X(t))dt + g(X(t)) dW(t),\;\; X(0)=X_0, \; 0 \le t \le T \ \ \ \ \ (*)
\]
Comments: 
1. 讀者需注意我們不可寫 $dW(t)/dt$ 因為 Browian motion 為 處處連續但處處不可微 (with probability 1)
2. 若 $g = 0$ 且 $X_0$ 為常數,則 SDE 退化成一般的 ODE;亦即
\[
dX(t) = f(X(t))dt, \;\; X(0)=X_0, \; 0 \le t \le T
\] (此時即可用 Euler method 求數值解。)
3. 關於 SDE 何時有解,讀者可參考BLOG 相關系列文章 :
[隨機分析] Uniqueness and Existence theorem for S.D.E. (1)- Uniqueness
[隨機分析] Uniqueness and Existence theorem for S.D.E. (2) - Picard Iteration for SDE
[隨機分析] Uniqueness and Existence theorem for S.D.E. (3) - An intermediate result (Upper bound) of Picard Iteration
[隨機分析] Uniqueness and Existence theorem for S.D.E. (4) - The Existence of Solution



那麼現在回歸主題,在此我們的目的是:
對 SDE 在關心的區間 $t \in [0,T]$ 之間進行數值求解 (透過 MATLAB )。

第一步首先將此區間 $[0,T]$ 離散化,亦即 對某些 $N$ 而言,定義 $\Delta t := T/N$ 且 對 $j=1,...,N$ ,定義 $\tau_j := j \Delta t$。 並且稱 數值近似解 $X_j :=X(\tau_j) $。

利用 Euler-Maruyama (EM) 法,我們可將前述 SDE 在離散化區間內 改寫如下\[{X_j} = {X_{j - 1}} + f({X_{j - 1}})\Delta t + g({X_{j - 1}})(W({\tau _j}) - W({\tau _{j - 1}})),\;\;j = 1,2,...,N
\]
Comment: 注意到上式若寫成積分形式即為
\[X({\tau _j}) = X({\tau _{j - 1}}) + \int_{{\tau _{j - 1}}}^{{\tau _j}} f (X(s))ds + \int_{{\tau _{j - 1}}}^{{\tau _j}} g (X(s))dW(s)
\]

Example
利用 Euler-Maruyama 法求解下列線性 SDE
\[
dX(t) = \mu X(t) dt + \sigma X(t) dW(t),\;\; X(0) = X_0
\]其中 $\mu, \sigma $ 為 常數 (亦即在此例我們選 $f(X(t)) = \mu X(t)$ 且 $g(X(t)) = \sigma X(t)$)。此例為 Geometric Brownian Motion (GBM) 模型,且此 SDE 有解析解如下:
\[X\left( t \right) = {X_0}{e^{\left( {\mu  - \frac{{{\sigma ^2}}}{2}} \right)t + \sigma W\left( t \right)}}\]以下我們將利用 MATLAB 實現 上述 解析解 與透過 EM 法 求數值解並比較其差異。


現在我們利用 MATLAB 實現 Euler-Maruyama 法:

首先計算 discretized Brownian path over $[0,1]$
定義 stepsize $\Delta t = R \delta t$ (一般泛取 $R=1$)且 由於 EM 法\[{X_j} = {X_{j - 1}} + f({X_{j - 1}})\Delta t + g({X_{j - 1}})(W({\tau _j}) - W({\tau _{j - 1}})),\;\;j = 1,2,...,L\]上式還需要計算 $W(\tau_j) - W(\tau_{j-1})$,故我們計算:
\[\begin{array}{l}
W({\tau _j}) - W({\tau _{j - 1}}) = W\left( {j\Delta t} \right) - W\left( {\left( {j - 1} \right)\Delta t} \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}&{}&{}&{}&{}
\end{array} = W\left( {jR\delta t} \right) - W\left( {\left( {j - 1} \right)R\delta t} \right)\\
\begin{array}{*{20}{c}}
{}&{}&{}&{}&{}&{}&{}&{}
\end{array} = \sum\limits_{k = jR - R + 1}^{jR} {d{W_k}}
\end{array}\]上式計算出現在下列程式碼 (第14行)


我們將 MATLAB 程式碼附上如下: (考慮 $\mu=2$, $\sigma =1$ 且 $X_0=1$。 $R=1$)
(line 10 : line 14: 產生 standard Brownian Motion)
(line 16: 解析解 for GBM SDE)
(line 19: line 25: EM-method)

執行結果如下:


上圖中 紅線為 EM 法結果藍線為 解析解結果。讀者可調整 $R$ 值來檢驗兩者差距。


ref: Desmond J. Higham, An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations, SIAM REVIEW Vol 43, No. 3, pp. 525-546, 2001.