1
12. Trigonometrická interpolace, rychlá Fourierova transformace, diskrétní Fourierova transformace Úloha trigonometrické interpolace je podobná jako u polynomiální interpolace. Budeme hledat trigonometrický polynom, tj. polynom tvaru
Pn (x) =
1 A0 + A1 cos x + B1 sin x + A2 cos 2x + B2 sin 2x + ... + An cos nx + Bn sin nx = 2
=
n X 1 (Aj cos jx + Bj sin jx), A0 + 2 j=1
který v daných uzlech nabývá stejných hodnot jako zadaná funkce. Trigonometrický polynom je funkce periodická s periodou 2π, budeme uvažovat zadanou funkci také periodickou se stejnou periodou. Budeme uvažovat sudý počet ekvidistantních uzlů, čili rozdělíme interval h0, 2πi na 2n intervalů stejné délky. Uzlové body označíme
ak = k ·
2π kπ = , 2n n
k = 0, 1, 2, ..., 2n − 1.
Lze ukázat, že v tomto případě má trigonometrický interpolační polynom tvar Pn (x) = 21 A0 + A1 cos 1x + B1 sin 1x + ... + An−1 cos(n − 1)x + Bn−1 sin(n − 1)x + 21 An cos nx =
=
n−1 X 1 1 (Aj cos jx + Bj sin jx) + An cos nx, A0 + 2 2 j=1
kde
Aj =
1 (f (a0 ) cos ja0 + f (a1 ) cos ja1 + ... + f (a2n−1 ) cos ja2n−1 ) = n 2n−1 1 X = f (ak ) cos jak , j = 0, 1, ..., n, n k=0
Bj =
1 (f (a0 ) sin ja0 + f (a1 ) sin ja1 + ... + f (a2n−1 ) sin ja2n−1 ) = n
=
2n−1 1 X f (ak ) sin jak , j = 0, 1, ..., n − 1. n k=0
1
2
P ř í k l a d 1: Stanovme trigonometrický interpolační polynom pro funkci, která je periodická s periπ odou 2π a v uzlových bodech má hodnoty: f (a0 ) = f (0) = 0, f (a1 ) = f ( ) = 1, f (a2 ) = f (π) = 2 3π 0, f (a3 ) = f ( ) = 1. 2 Interval h0, 2πi je rozdělen na čtyři intervaly, tedy n = 2. Vypočteme koeficienty A0 = A1 = A2 = B1 =
1 2 1 2 1 2 1 2
3 P
f (ak ) cos 0ak = 12 (0 + 1 + 0 + 1) = 1,
3 P
f (ak ) cos 1ak = 12 (0 + 1 · cos π2 + 0 + 1 · cos 3π 2 ) = 0,
3 P
2·3π f (ak ) cos 2ak = 12 (0 + 1 · cos 2π 2 + 0 + 1 · cos 2 ) = −1,
3 P
f (ak ) sin 1ak = 12 (0 + 1 · sin π2 + 0 + 1 · sin 3π 2 )=0
k=0
k=0
k=0
k=0
a dosadíme do polynomu P2 (x) =
1 1 1 1 + 0 + 0 + · (−1) · cos 2x = − cos 2x. 2 2 2 2
Rychlá Fourierova transformace Rychlá Fourierova transformace je název pro algoritmus, který slouží pro rychlý výpočet součtů tvaru N −1 X
e
2πiqs N
fs , q = 0, ..., N − 1,
s=0
kde fs jsou komplexní čísla. Výraz eiφ je definován následovně: eiφ = cos φ + i sin φ. Neboť funkce kosinus je sudá a sinus lichá, platí: e−iφ = cos φ − i sin φ, Sečtením a odečtením uvedených rovností získáme eiφ + e−iφ = 2 cos φ
odtud
cos φ =
eiφ − e−iφ = 2 sin φ
odtud
sin φ =
eiφ +e−iφ . 2 eiφ −e−iφ . 2
Koeficienty trigonometrické interpolace lze tedy vypočítat podle vztahů Aj = Cj + C−j , Bj = Cj − C−j , kde Cj =
2n−1 1 X f (ak )eijak , j = 0, ±1, ..., ±n. 2n k=0
Nyní budeme počítat součty N −1 X
e
2πiqs N
fs , q = 0, ..., N − 1.
s=0
2πiq
Uvedený součet lze zapsat ve tvaru pN −1 (z) = f0 + f1 z + ... + fN −1 z N −1 , kde z = e N , q = 0, 1, ..., N − 1. Jde o polynom stupně N − 1, jehož hodnoty máme spočítat v N − 1 bodech. Předpokládejme, že N = 2k pro jisté reálné číslo k. Pak lze psát 2
3
pN −1 (z) = zp N −1 (z 2 ) + r N −1 (z 2 ), 2
kde p N −1 , r N −1 jsou polynomy stupně 2
2
N 2
2
− 1,
N 2
p N −1 (u) = 2
−1 X
N 2
j
f2j+1 u ,
r N −1 (u) = 2
j=0
−1 X
f2j uj .
j=0
Nyní budeme počítat hodnoty dvou polynomů stupně 2πiq 2 2πi(q− N2 ) 2 N q = 2 , ..., N − 1 je e N = e N .
N 2
− 1 v stupně
N 2
bodech, neboť pro
Dále podobně vyjádříme polynomy p N −1 , r N −1 jako součty polynomů stupně N4 − 1, a tak dále 2 2 až dojdeme k polynomům stupně 1, jejichž hodnoty budeme počítat v bodech ±1. Potřebný počet operací nepřevýší číslo 2N log2 N , zatímco při přímém výpočtu koeficientů je třeba řádově N 2 operací. P ř í k l a d 2: Stanovme trigonometrický interpolační polynom pro funkci, která je periodická s periπ odou 2π a v uzlových bodech má hodnoty: f (a0 ) = f (0) = 0, f (a1 ) = f ( ) = 1, f (a2 ) = f (π) = 2 3π 0, f (a3 ) = f ( ) = 1. Použijeme rychlou Fourierovu transformaci. V tomto případě to není rychlejší 2 postup, protože se jedná o polynomy nízkého stupně. Budeme počítat koeficienty Cj = Označme Dj =
3 P
1 4
3 P
f (ak )eijak , j = 0, ±1, ±2.
k=0
f (ak )eijak , j = 0, ±1, ±2, fk = f (ak ).
k=0
Konkrétně: 0·iπ 0 0·iπ 1 0·iπ 2 0·iπ 3 D0 = f 0 e 2 + f1 e 2 + f2 e 2 + f3 e 2 = f0 z 0 + f1 z 1 + f2 z 2 + f3 z 3 = f0 + f2 z 2 + 0·iπ 2 0·iπ z(f1 + f3 z 2 ) = r1 (u) + zp1 (u), pro u = e 2 = 1, z = e 2 = 1, 1·iπ 1 1·iπ 2 1·iπ 3 1·iπ 0 + f1 e 2 + f2 e 2 + f3 e 2 = f0 z 0 + f1 z 1 + f2 z 2 + f3 z 3 = f0 + f2 z 2 + D1 = f 0 e 2 1·iπ 2 1·iπ z(f1 + f3 z 2 ) = r1 (u) + zp1 (u), pro u = e 2 = −1, z = e 2 = i,
−1·iπ 0 −1·iπ 1 −1·iπ 2 −1·iπ 3 D−1 = f0 e 2 + f1 e 2 + f2 e 2 + f3 e 2 = f0 z 0 + f1 z 1 + f2 z 2 + f3 z 3 = −1·iπ 2 −1·iπ f0 + f2 z 2 + z(f1 + f3 z 2 ) = r1 (u) + zp1 (u), pro u = e 2 = −1, z = e 2 = −i,
2·iπ 1 2·iπ 2 2·iπ 3 2·iπ 0 + f1 e 2 + f2 e 2 + f3 e 2 = f0 z 0 + f1 z 1 + f2 z 2 + f3 z 3 = f0 + f2 z 2 + D2 = f 0 e 2 2·iπ 2 2·iπ z(f1 + f3 z 2 ) = r1 (u) + zp1 (u), pro u = e 2 = 1, z = e 2 = −1,
−2·iπ 0 −2·iπ 1 −2·iπ 2 −2·iπ 3 D−2 = f0 e 2 + f1 e 2 + f2 e 2 + f3 e 2 = f0 z 0 + f1 z 1 + f2 z 2 + f3 z 3 = −2·iπ 2 −2·iπ f0 + f2 z 2 + z(f1 + f3 z 2 ) = r1 (u) + zp1 (u), pro u = e 2 = 1, z = e 2 = −1. Vidíme, že hodnoty u se začínají opakovat.
Určíme polynomy r1 (u) = 0 + 0u = 0, p1 (u) = 1 + 1u, p1 (1) = 2, p1 (−1) = 0, . Po dosazení dostaneme: 3
4
D0 = 0 + 1 · p1 (1) = 2, D1 = 0 + i · p1 (−1) = 0, D−1 = 0 − i · p1 (−1)) = 0, D2 = 0 − 1 · p1 (1) = −2, D−2 = 0 − 1 · p1 (1)) = −2. Tedy C0 = 14 , D0 = 21 , podobně C±1 = 0, C±2 = − 12 . Nakonec dosadíme do vzorců Aj = Cj + C−j , Bj = Cj − C−j , dostaneme A0 =
1 , A1 = 0, A2 = −1, B1 = 0. 2
Trigonometrický polynom získáme dosazením jako v předchozím příkladě. Diskrétní Fourierova transformace Konečná posloupnost F(a) = (b0 , b1 , ..., bN −1 ), kde bq =
N −1 X
e−
2πiqs N
as , q = 0, ..., N − 1
s=0
se nazývá diskrétní Fourierova transformace posloupnosti a = (a0 , a1 , ..., aN −1 ). Naopak posloupnost F −1 (b) = (a0 , a1 , ..., aN −1 ) se nazývá zpětná (inverzní) diskrétní Fourierova transformace posloupnosti b = (b0 , b1 , ..., bN −1 ) a je dána vztahy as =
N −1 1 X 2πiqs e N bq , s = 0, ..., N − 1. N q=0
Z názvu je patrné, že diskrétní Fourierovou transformací zpětné diskrétní Fourierovy tranformace získáme opět původní posloupnost, tedy F −1 (F(a)) = a, Navíc platí: F −1 (c) =
1 (F(rev(c))), N
F(F −1 (b)) = b.
F(c) = N (F −1 (rev(c))),
kde rev(c) = (c0 , cn−1 , cn−2 , ..., c2 , c1 ). K výpočtu diskrétní Fourierovy transformace slouží algoritmus rychlé Fourierovy transformace. P ř í k l a d 3: Určeme diskrétní Fourierovu transformaci posloupnosti a = (0, 1, 0, 1). 2π
s=0
,
2π
π
Označme z = (z0 , z1 , z2 , z3 ), kde zi = e− N = e− 4 = e− 2 . Tedy z = (1, −i, −1, i). Pak 3 3 X X 2πiqs e− 4 as = bq = (zsq )s as , q = 0, ..., 3 s=0
dostaneme b0 = (z00 )0 a0 + (z10 )0 a1 + (z20 )0 a2 + (z30 )0 a3 = 0 + 1 + 0 + 1 = 2, 4
5
b1 = (z01 )0 a0 + (z11 )0 a1 + (z21 )0 a2 + (z31 )0 a3 = 0 − i + 0 + i = 0, b2 = (z02 )0 a0 + (z12 )0 a1 + (z22 )0 a2 + (z32 )0 a3 = 0 − 1 + 0 − 1 = −2, b3 = (z03 )0 a0 + (z13 )0 a1 + (z23 )0 a2 + (z33 )0 a3 = 0 + i + 0 − i = 0. Stačí si uvědomit, že první tři členy hledané posloupnosti jsou čísla D0 , D−1 , D−2 , které jsme vypočetli v Příkladu 2. Výpočet čtvrtého členu by probíhal podobně, pro úplnost uvádíme, že D−3 = 0. Pro vektory (konečné posloupnosti) a, b definujeme konvoluci vektorů a, b jako vektor c = a ∗ b,, kde N −1 X aj−k bk , cj = q=0
kde pro záporné indexy j − k definujeme j − k = n + j − k. P ř í k l a d 4: (1, 2, 0, 1, 2) ∗ (0, −1, 0, 1, 1) = (0, 0, 1, 3, 2). Oceňování opcí pomocí diskrétní Fourierovy transformace Uvažujme binomický model oceňování opcí, kde C(j) je vektor cen opce v časej = 0, 1, ..., N , q = (qu , qd , 0, ..., 0) je vektor (risk-neutral) pravděpodobností doplněný nulami na dalších místech, 1 r ) 12 , kde r je (risk-free) úroková míra. Pak platí: rev(q) = (qu , 0, ..., 0, qd ) , Rf = (1 + 100 C(N − 1) = C(N ) ∗ rev(q)/Rf , C(N − 2) = C(N − 1) ∗ rev(q)/Rf = C(N ) ∗ rev(q) ∗ rev(q)/Rf2 , ...... (N −j)krát
z }| { C(j) = C(N ) ∗ rev(q) ∗ rev(q) /Rfj . V C(j) uvažujeme pouze prvních j + 1 hodnot. Tedy N krát
z }| { C(j) = C(N ) ∗ rev(q) ∗ rev(q) /Rfj . Počáteční hodnota opce (v čase 0) je rovna první složce tohoto vektoru. P ř í k l a d 5: V binomickém modelu jsme zjistili, že qu = 0.43523, qd = 0.56477, r = 4%, tedy Rf = 1.0033. 0.43523Cu + 0.56477Cd qu Cu + qd Cd = . Platí: hodnota value(C) = Rf 1.0033 Hodnota opce je dána tabulkou (splatnost za 3 měsíce): 5
6
number of low returns
C(0)
C(1)
C(2)
C(3)
0
81.36
162.66
317.54
599.64
19.19
44.25
102.00
0
0
1 2 3
0
Hodnoty počítáme postupně, např. 317, 54 =
0, 43523 · 599, 64 + 0, 56477 · 102 . 1, 0033
Sloupec C(2) lze též vypočítat jako konvoluci 0, 56477 0, 43523 ; 0; 0; )= C(3) ∗ rev(q/Rf ) = (599, 6; 102; 0; 0) ∗ ( 1, 0033 1, 0033 = (599, 6; 102; 0; 0) ∗ (0, 43523; 0; 0; 0, 56477) = (317, 5400; 44, 2474; 0; 337, 5448). Bereme v úvahu první tři složky vypočteného vektoru. Diskrétní Fourierova transformace má užitečnou vlastnost: F(a ∗ b) = F(a)F(b), kde součin na pravé straně je vektor, jehož složky jsou součiny složek daných vektorů. Pak konvoluci lze vypočítat následovně: a ∗ b = F −1 (F(a)F(b)). Tudíž (nearbitrážní) cena opce je rovna v čase 0 je rovna první složce vektoru 1 C(j) = F −1 F(C(N ))F( rev(q))N . Rf MATLAB - rychlá Fourierova transformace Existuje příkaz y = f f t(x), resp. y = if f t(x), který určí diskrétní Fourierovu, resp. inverzní F. transformaci posloupnosti (vektoru) x, která je počítána pomocí rychlé Fourierovy transformace. Pokud bude vstup ve tvaru matice X, bude výsledkem diskrétní Fourierova transformace každého jejího sloupce.
6
