Εισαγωγή στις Σ.Δ.Ε.

Unprotect[ToString];
ToString[expr: _[__], StandardForm] := With[{view = MMAView[expr]}, ExportString[
    StringReplace[
        (view // ToBoxes) /. {RowBox->RowBoxFlatten} // ToString
    , {"\[NoBreak]"->""}]
, "String"]];
Protect[ToString];
Clear["Global`*"]

Θεμελίωση της έννοιας

Μια απλή Διαφορική Εξίσωση είναι η $x'=b(x,t)$ ή, κάνωντας χρήση διαφορικών, $dx=b(x,t)dt$. Αυτό μάς λέει με ποιον τρόπο μεταβάλλεται το $x$ μετά από κάποιο μικρό χρονικό διάστημα $dt$.

Στην πραγματική ζωή, όμως, τίποτα δεν είναι απόλυτα ντετερμινιστικό. Τουλάχιστον βάσει του επιπέδου γνώσης που έχουμε ή που μπορούμε να διαχειριστούμε. Έτσι, μια νότα τυχαιότητας στην παραπάνω Δ.Ε. θα έδινε ένα άρωμα ρεαλισμού εντονότερο από πριν. Μια ιδέα είναι να αυξομοιώνεται το αποτέλεσμα βάσει ενός τυχαίου μικρού άλματος $dW_t$ διεσταλμένο βάσει μιας συνάρτησης $\sigma(X_t, t)$. Έτσι έχουμε την παρακάτω Στοχαστική Διαφορική Εξίσωση και λέμε ότι η $X_t$ είναι μία διαδικασία Itô:

\[ dX_t = b(X_t, t)dt+\sigma(X_t, t)dW_t \]

Στον παραπάνω τύπο το $W_t$ είναι μια διαδικασία Wiener (βλ. [2] σελ. 284 και [5] σελ. 290). Δηλαδή έχει τις κάτωθι ιδιότητες (βλ [3] σελ. 608):

Η διαδικασία Itô μπορεί να εκφραστεί και ως:

\[ X_t=x_0+\int_0^t b(X_s, s)ds+\int_0^t \sigma(X_s, s)dW_s, \]

όπου το $x_0$ είναι η γνωστή αρχική τιμή της διαδικασίας (την αρχική τιμή τη γνωρίζουμε, δεν είναι τυχαία) και το $\int_0^t \sigma(X_s, s)dW_s$ είναι ένα στοχαστικό ολοκλήρωμα.

Αν μία φορά τα ολοκληρώματα οι Διαφορικές Εξισώσεις ήταν ζόρικα πεδία και απαιτούσαν ειδικά τεχνάσματα, τα στοχαστικά ολοκληρώματα και οι Στοχαστικές Διαφορικές Εξισώσεις είναι για να απελπίζεσαι. Έτσι, δεν θα δώσουμε εδώ μια γενική λύση της παραπάνω Σ.Δ.Ε., αλλά θα λύσουμε μόνο κάποιες υποπεριπτώσεις της, και στη γενικότερη περίπτωση θα αρκεστούμε στη γραφική της επίλυση χρησιμοποιώντας μια διακριτοποιημένη εκδοχή της διαδικασίας Itô (βλ. [5] σελ. 291):

\[ \Delta X=b(X_t,t)\Delta t+ \sigma(X_t,t)Z\sqrt{\Delta t}, \]

όπου $\Delta X=X_{t+\Delta t}-X_t$ και $Z\sim N(0,1)$. Ενδεικτικά παρουσιάζουμε τη γραφική αναπαράσταση της Σ.Δ.Ε.:

\[ dX_t = (5X_t+1.5t)dt+(1.1x+1-e^{-t})dW_t \]

EulerMaruyamaPlot[x0_, b_, sigma_, dt_, T_, n_, title_] := 
 Module[{times, nSteps, data},
  
  (* Δημιουργία του διακριτού χρονικού άξονα *)
  times = Range[0, T, dt];
  nSteps = Length[times] - 1;
  
  (* Αν το T δεν είναι ακριβές πολλαπλάσιο του dt, το διορθώνουμε *)
  If[Abs[times[[-1]] - T] > 10^-10, 
   times = Range[0, nSteps*dt, dt];
  ];
  
  (* Προσομοίωση n ανεξάρτητων τροχιών *)
  data = Table[
    Module[{path},
      path = FoldList[
        #1 + b[#1, #2]*dt + sigma[#1, #2]*Sqrt[dt]*
          RandomVariate[NormalDistribution[0, 1]] &,
        x0,
        Most[times]
      ];
      Transpose[{times, path}]
    ],
    {n}
  ];
  
  (* Γραφική παράσταση με τις βελτιώσεις *)
  ListLinePlot[data,
   PlotRange -> All,
   AxesLabel -> {"t", "X(t)"},
   PlotLabel -> 
    title <> "\n(x0 = " <> ToString[x0] <> ", dt = " <> ToString[dt] <> 
     ", T = " <> ToString[T] <> ", n = " <> ToString[n] <> ")",
   Epilog -> {
     (* Κόκκινη κουκκίδα στο σημείο (0, x0) *)
     PointSize[Large],
     Red,
     Point[{0, x0}]
   },
   ImageSize -> Large,
   PlotTheme -> "Scientific"
  ]
 ]
b1[x_, t_] := 5 x+1.5*t;
sigma1[x_, t_] := 1.1 x+1-Exp[-t];

EulerMaruyamaPlot[100, b1, sigma1, 0.01, 1, 10, "Τροχιές Σ.Δ.Ε."]
Plot

Το λήμμα Itô

Σημαντικό για την επίλυση κάποιων Σ.Δ.Ε., αλλά και για τον υπολογισμό διαφόρων χαρακτηριστικών της $X_t$ είναι το λήμα Itô (βλ. [1] σελ. 383 και [5] σελ. 292), το οποίο δίνει το διαφορικό μιας συνάρτησης ($G$) των $X_t$ και $t$. Έτσι για την πραγματική συνάρτηση $G(x,t)$ με συνεχείς τις μερικές παραγώγους της ($\partial_t G$, $\partial_x G$ και $\partial_{xx}G$) έχουμε:

\[\begin{align*} dG(X_t,t)=&\left(\dfrac{\partial G(X_t,t)}{\partial t}+b(X_t,t)\dfrac{\partial G(X_t,t)}{\partial x}+\dfrac{\sigma^2(X_t,t)}{2}\dfrac{\partial^2 G(X_t,t)}{\partial x^2}\right)dt\\ &+\sigma(X_t,t)\dfrac{\partial G(X_t,t)}{\partial x}dW_t \end{align*}\]

Στις ακόλουθες υποενότητες θα χρησιμοποιήσουμε το λήμμα Itô για να επανεξετάσουμε κάποιους κανόνες του γνωστού μας Διαφορικού Λογισμού.

Παράγωγος γινομένου

Εχουμε ότι:

\[ dG(X_t)=\left(b(X_t,t)G'(X_t)+\dfrac{\sigma^2(X_t,t)}{2}G''(X_t)\right)dt+\sigma(X_t,t)G'(X_t)dW_t \]

και άρα ότι:

\[\begin{align*} d\left(G(X_t)H(t)\right)=&\left(G(X_t)H'(t)+b(X_t,t)G'(X_t)H(t)+\dfrac{\sigma^2(X_t,t)}{2}G''(X_t)H(t)\right)dt\\ &+\sigma(X_t,t)G'(X_t)H(t)dW_t\\ =& \ G(X_t)H'(t)dt+\left(b(X_t,t)G'(X_t)H(t)+\dfrac{\sigma^2(X_t,t)}{2}G''(X_t)H(t)\right)dt\\ &+\sigma(X_t,t)G'(X_t)H(t)dW_t\\ =& \ G(X_t)dH(t)+H(t)dG(X_t) \end{align*}\]

Η τελευταία σχέση είναι η γνωστή μας παραγώγιση γινομένου και θα μάς χρησιμεύσει στην επίλυση της διαδικασίας Ornstein-Uhlenbeck. Ο αναγνώστης αξίζει να προσέξει ότι για την «επιβίωση» του γνωστού κανόνα παραγώγισης γινομένου έπρεπε ο ένας όρος (το $H(t)$) να είναι ντετερμινιστικός.

Παράγωγος σύνθεσης

Όπως πριν, έχουμε ότι:

\[ dG(X_t)=\left(b(X_t,t)G'(X_t)+\dfrac{\sigma^2(X_t,t)}{2}G''(X_t)\right)dt+\sigma(X_t,t)G'(X_t)dW_t \]

Άρα:

\[\begin{align*} d\left(H\left(G(X_t)\right)\right)=& \ b(X_t,t)H'\left(G(X_t)\right)G'(X_t)dt\\&+\left(\dfrac{\sigma^2(X_t,t)}{2}\left(H''\left(G(X_t)\right)G'(X_t)^2+H'\left(G(X_t)\right)G''(X_t)\right)\right)dt\\ &+\sigma(X_t,t)H'\left(G(X_t)\right)G'(X_t)dW_t\\ =& \ H'\left(G(X_t)\right)\left(b(X_t,t)G'(X_t)+\dfrac{\sigma^2(X_t,t)}{2}G''(X_t)\right)dt\\ &+\sigma(X_t,t)H'\left(G(X_t)\right)G'(X_t)dW_t\\ &+\dfrac{\sigma^2(X_t,t)}{2}H''\left(G(X_t)\right)G'(X_t)^2 dt\\ =& \ H'\left(G(X_t)\right)dG(X_t)+\dfrac{\sigma^2(X_t,t)}{2}H''\left(G(X_t)\right)G'(X_t)^2 dt \end{align*}\]

Εδώ βλέπουμε πως παύει πλέον να ισχύει ο κανόνας $d\left(H\left(G(X_t)\right)\right)=H'\left(G(X_t)\right)dG(X_t)$, εκτός κι αν η $H$ είναι γραμμική, οπότε $H''(G(X_t))=0$.

Μέση τιμή και διασπορά

Σε κάθε στοχαστική διαδικασία επιχειρούμε να προσδιοίσουμε κάποια στατιστικά μεγέθη της, τα οποία θα συνοψίζαν την πιθανοκρατική συμπεριφορά της. Έτσι θα πράξουμε κι εδώ, δίνοντας γενικές οδηγίες περί του πώς μπορούμε να βρούμε:

μιας διαδικασίας Itô.

Δεδομένου ότι η $X_t$ και το απειροστό $dW_t$ είναι ανεξάρτητα, όπως και τα $dW_t$ μεταξύ τους, έχουμε από την $dX_t = b(X_t, t)dt+\sigma(X_t, t)dW_t$:

\[\begin{align*} &\mathbb{E}(dX_t) = \mathbb{E}\left(b(X_t, t)dt+\sigma(X_t, t)dW_t\right)\\ \Leftrightarrow \ & d\mathbb{E}X_t=\mathbb{E}b(X_t, t)dt+\mathbb{E}\left(\sigma(X_t, t)dW_t\right)\\ \Leftrightarrow \ & d\mathbb{E}X_t=\mathbb{E}b(X_t, t)dt+\left(\mathbb{E}\sigma(X_t, t)\right)\left(d\mathbb{E}W_t\right)\\ \Leftrightarrow \ & d\mathbb{E}X_t=\mathbb{E}b(X_t, t)dt \end{align*}\]

όπου η τελευταία ισότητα τεκμέρεται από το γεγονός $\mathbb{E}W_t=0$. Η σχέση $d\mathbb{E}X_t=\mathbb{E}b(X_t, t)dt$, που φτιάξαμε, δεν εμπεριέχει πουθενά στοιχεία Στοχαστικού Λογισμού. Είναι μια απλή Διαφορική Εξίσωση με άγνωστη συνάρτηση την $m(t)=\mathbb{E}X_t$. Ο αναγνώστης αξίζει να παρατηρήσει ότι πρόκειται για τη Διαφορική Εξίσωση που σχηματίζεται αν αγνωήσουμε τον όρο που προσδίδει πιθανοκρατική συμπεριφορά στην $X_t$, τον $\sigma(X_t, t)dW_t$ και κρατήσουμε μόνο το ντετερμινιστικό μέρος.

Για τον υπολογισμό της $\operatorname{Var}(X_t)$, λόγω της:

\[ \operatorname{Var}(X_t)=\mathbb{E}X_t^2-(EX_t)^2 \]

και δεδομένου ότι έχουμε βρει την $\mathbb{E}X_t$, δεν έχουμε, παρά να βρούμε την $\mathbb{E}X_t^2$. Έτσι, λόγω του λήμματος Itô (όπου $G(x,t)=x^2$) έχουμε:

\[ dX_t^2=\left(2b(X_t,t)X_t+\sigma^2(X_t,t)\right)dt+2\sigma(X_t,t)X_tdW_t \]

Έτσι, όπως πριν:

\[\begin{align*} &\mathbb{E}(dX_t^2) = \mathbb{E}\left(\left(2b(X_t,t)X_t+\sigma^2(X_t,t)\right)dt+2\sigma(X_t,t)X_tdW_t\right)\\ \Leftrightarrow & \ d\mathbb{E}X_t^2=\mathbb{E}\left(2b(X_t,t)X_t+\sigma^2(X_t,t)\right)dt+\mathbb{E}\left(2\sigma(X_t,t)X_tdW_t\right)\\ \Leftrightarrow & \ d\mathbb{E}X_t^2=\mathbb{E}\left(2b(X_t,t)X_t+\sigma^2(X_t,t)\right)dt+\mathbb{E}\left(2\sigma(X_t,t)X_t\right)\mathbb{E}\left(dW_t\right)\\ \Leftrightarrow & \ d\mathbb{E}X_t^2=\left(\mathbb{E}\left(2b(X_t,t)X_t\right)+\mathbb{E}\left(\sigma^2(X_t,t)\right)\right)dt \end{align*}\]

Η $d\mathbb{E}X_t^2=\left(\mathbb{E}\left(2b(X_t,t)X_t\right)+\mathbb{E}\left(\sigma^2(X_t,t)\right)\right)dt$ μια φυσιολογική Διαφορική Εξίσωση που μάς φανερώνει τη συνάρτηση $sq(t)=\mathbb{E}X_t^2$, ήπερ είναι το τελευταίο ληθαράκι, για να βρούμε την $\operatorname{Var}(X_t)$.

Ας δούμε τα παραπάνω με ένα παράδειγμα μιας Σ.Δ.Ε. Σε πρώτη φάση ας ρίξουμε μια ματιά στη συνήθη Διαφορική Εξίσωση:

\[ x'=\alpha\cdot(\mu-x), \ \alpha>0. \]

Η φυσική της σημασία είναι προφανής: Η μεταβλητή $x$ έχει την τάση να επιστρέφει στο $\mu$, διότι:

Με διαφορικά η εν λόγω εξίσωση γράφεται $dx=\alpha(\mu-x)dt$. Θέλοντας να κρατήσουμε την προαναφερθείσα συλλογιστική επαναφοράς, αλλά να τη διαταράξουμε με λίγη τυχαιότητα, τότε έχουμε την εξίσωση Langevin, η λύση της οποίας καλέιται διαδικασία Ornstein-Uhlenbeck (βλ. [5] σελ. 487-489):

\[ dX_t=\theta\cdot(\mu-X_t)dt+\sigma dW_t, \ \alpha,\sigma>0. \]

Παρουσιάζουμε κάποιες εκδοχές της για επαλήθευση της συμπεριφοράς που οικοδομήσαμε να έχει:

bOU[x_, t_] := 2(3-x);
sigmaOU[x_, t_] := 0.2;

EulerMaruyamaPlot[1, bOU, sigmaOU, 0.01, 3, 10,"Τροχιές Ornstein-Uhlenbeck με θ=2, μ=3 και σ=0.2"]
EulerMaruyamaPlot[5, bOU, sigmaOU, 0.01, 3, 10,"Τροχιές Ornstein-Uhlenbeck με θ=2, μ=3 και σ=0.2"]
Plot
Plot

Έχουμε, όπως πριν:

\[ \mathbb{E}\left(dX_t\right)=\mathbb{E}\left(\theta\cdot(\mu-X_t)dt+\sigma dW_t\right)\Leftrightarrow d\mathbb{E}X_t=\theta(\mu-EX_t)dt \]

Θέτουμε $m(t)=\mathbb{E}X_t$, άρα έχουμε τη Διαφορική Εξίσωση:

\[ m'=\theta (\mu-m) \]

με $m(0)=\mathbb{E}X_0=x_0$, την οποία και λύνουμε.

solEou = DSolve[{m'[t]==θ*(μ-m[t]),m[0]==x0},m[t],t]//Expand//Flatten;
solMou = m[t]/.solEou
Plot

Βρήκαμε ότι:

\[ m(t)=\mu+e^{-\theta t}\left(x_0-\mu\right) \]

Ένα πρώτο συμπέρασμα είναι ότι καθώς $t\to+\infty$ είναι $\mathbb{E}X_t\to \mu$. Το αν και η ίδια η $X_t$ συγκλίνει στη $\mu$ ή αν η τυχαιότητα την στέλνει σε όλα τα σημεία του ορίζοντα, αυτό θα μάς το απαντήσει η $\operatorname{Var}(X_t)$.

Για τη διακύμανση θα πάμε να υπολογίσουμε το $dX_t^2$, όπως είπαμε και στις γενικές οδηγίες. Λόγω του λήμματος Itô (όπου $G(x,t)=x^2$) έχουμε:

\[ dX_t^2=\left(2\theta (\mu-X_t)X_t+\sigma^2\right)dt+2\sigma X_tdW_t. \]

«Κοπανάμε» $\mathbb{E}$ και στις δύο μεριές της παραπάνω εξίσωσης, οπότε έχουμε:

\[\begin{align*} &\mathbb{E}(dX_t^2) = \mathbb{E}\left(\left(2\theta (\mu-X_t)X_t+\sigma^2\right)dt+2\sigma X_tdW_t\right)\\ \Leftrightarrow \ & d\mathbb{E}X_t^2=\mathbb{E}\left(2\theta (\mu-X_t)X_t+\sigma^2\right)dt+\mathbb{E}\left(2\sigma X_tdW_t\right)\\ \Leftrightarrow \ & d\mathbb{E}X_t^2=\left(2\theta\mu\mathbb{E}X_t-2\theta\mathbb{E}X_t^2+\sigma^2\right)dt\\ \Leftrightarrow \ & d\mathbb{E}X_t^2=\left(2\theta\mu\left(\mu+e^{-\theta t}\left(x_0-\mu\right)\right)-2\theta\mathbb{E}X_t^2+\sigma^2\right)dt \end{align*}\]

Θέτοντας $sq(t)=\mathbb{E}X_t^2$ έχουμε τη συνήθη Διαφορική Εξίσωση:

\[ sq'(t)=2\theta\mu\left(\mu+e^{-\theta t}\left(x_0-\mu\right)\right)-2\theta sq(t)+\sigma^2 \]

την οποία και θα λύσουμε, θεωρώντας πάλι $sq(0)=\mathbb{E}X_0^2=x_0^2$.

solVarOU = DSolve[{sq'[t]==2θ*μ*(μ+E^(-θ*t)(x0-μ))-2θ sq[t]+σ^2,sq[0]==x0^2},sq[t],t][[1]]//Expand//Simplify;
solSQou = sq[t]/.solVarOU
Plot

Λύνοντας την παραπάνω Διαφορική Εξίσωση βρήκαμε:

\[ sq(t)=\frac{e^{-2\theta t}\left(\left(e^{\theta t}-1\right)\left(2\theta\mu^2\left(e^{\theta t}-1\right)+\sigma^2\left(e^{\theta t}+1\right)\right)+4\theta\mu x_0\left(e^{\theta t}-1\right)+2\theta x_0^2\right)}{2\theta} \]

Άρα έχουμε:

\[\begin{align*} \operatorname{Var}X_t&=\mathbb{E}X_t^2-(\mathbb{E}X_t)^2= \end{align*}\]
solSQou-solMou^2//FullSimplify
Plot

Βλέπουμε ότι η $\operatorname{Var}(X_t)$:

Την φραγή της τυχαιότητας τη διαπιστώνουμε και από τη γραφική αναπαράσταση κάποιων Ornstein-Uhlenbeck, όπου είναι εμφανές ότι η διασπορά των τιμών αυξάνεται ναι μεν, αλλά όχι πέρα από κάποιο σημείο. Κάνοντας έτσι κάθε υλοποίηση της Ornstein-Uhlenbeck να γυροφέρνει την καμπύλη της μέσης τιμής χωρίς μεγάλες διαφοροποιήσεις.

bOU2[x_, t_] := 2(5-x);
sigmaOU2[x_, t_] := 0.7;

EulerMaruyamaPlot[1, bOU2, sigmaOU2, 0.01, 3, 10,"Τροχιές Ornstein-Uhlenbeck με θ=2, μ=5 και σ=0.7"]
Plot

Η εξίσωση Fokker-Plank για τη σ.π.π.

Δεδομένης της μη-επιλυσιμότητας της Σ.Δ.Ε. πολλάκις αρκούμαστε στο να βρούμε την κατανομή της $X_t$ κι όχι την $X_t$ την ίδια. Η συνάρτηση πυκνότητας πιθανότητας της $X_t$ (έστω $p(x,t)$) είναι λύση της κάτωθι Μερικής Διαφορικής Εξίσωσης (βλ. [1] σελ. 380 και 382), όπου ο αναγνώστης θα την βρει στη βιβλιογραφία είτε ως εξίσωση Fokker-Plank είτε ως εξίσωση Kolmogorov ορθής φοράς:

\[ \dfrac{\partial p(x,t)}{\partial t}=-\dfrac{\partial\left(b(x,t)p(x,t)\right)}{\partial x}+\dfrac{1}{2}\dfrac{\partial^2\left(\sigma^2(x,t)p(x,t)\right)}{\partial x^2} \]

με:

Θα λύσουμε την εν λόγω Μερική Διαφορική Εξίσωση για την περίπτωση της Ornstein-Uhlenbeck.

⚠️ Warning
⚠️ ΥΠΟ ΚΑΤΑΣΚΕΥΗ ⚠️

Διάφορες λύσεις

Στην ενότητα αυτή θα παρουσιάσουμε διάφορα στοχαστικά μοντέλα κοινωνικών ή οικολογικών φαινομένων που να θεμελιώνονται με τη βοήθεια Σ.Δ.Ε. και θα προσπαθήσουμε να προσδιορίσουμε την ίδια τη $X_t$ κι όχι μόνο κάποιες στατιστικές μετρήσεις της, όπως η μέση τιμή και η διασπορά της.

Αντίστοιχα με τη στοχαστική παραγώγιση, έχουμε, όπως προαναφέραμε, και τη στοχαστική ολοκλήρωση. Δεν θα επεκταθούμε εδώ σε ορισμούς και εξηγήσεις, θα αναφέρουμε μόνο ότι ισχύουν οι βασικές ιδιότητες των ολοκληρωμάτων που ξέρουμε (προ Θ.Θ.Ο.Λ.), αλλά και η:

\[ \int_{t_1}^{t_2}dG(X_t,t)=G(X_{t_2},t_2)-G(X_{t_1},t_1) \]

Σημαντική είναι και η κάτωθι ιδιότητα (βλ. [5] σελ. 308-309 και [1] σελ. 384):

\[ \int_{t_1}^{t_2}W_sdW_s=\dfrac{W_{t_2}^2-W_{t_1}^2}{2}-\dfrac{t_2-t_1}{2} \]

Ο αναγνώστης αξίζει να παρατηρήσει τη διαφορά που έχει αυτή έναντι της γνωστής ντετερμινιστικής ολοκλήρωσης:

\[ \int_{t_1}^{t_2}W_t dW_t=\dfrac{W_{t_2}^2-W_{t_1}^2}{2}. \]

Γεωμετρική κίνηση Brown

Σχεδόν από το λύκειο γνωρίζει κανείς το εκθετικό πληθυσμιακό μοντέλο, σύμφωνα με το οποίο η αύξηση ενός πληθυσμού είναι ανάλογη του μεγέθους του. Ένα απόλυτα λογικό μοντέλο, καθόσον, αν ο καθένας από εμάς τους $x$ Έλληνες κάνει από $2$ παιδιά, τότε η αύξηση του πληθυσμού μας (τα συνολικά παιδιά, δηλαδή) θα είναι $\Delta x=2x$. Σε μαθηματικό συμβολισμό έχουμε:

\[ dx=rxdt. \]

Ας βάλουμε λίγο θόρυβο στο μοντέλο αυτό. Συγκεκριμένα, ας θεωρήσουμε ότι το εύρος των τυχαίων κινήσεων είναι ανάλογο του ληθυσμού, μια εύλογη υπόθεση, αν θεωρήσουμε ότι οι τυχαίες διακυμάνσεις οφείλονται σε κάποιους τυχαίους παράγοντες που αφορούν το κάθε ένα άτομο. Έτσι έχουμε την κάτωθι Σ.Δ.Ε. (βλ. [1] σελ. 384-385), η οποία ακούς στο όνομα γεωμετρική κίνηση Brown.

\[ dX_t=rX_tdt+cX_tdW_t. \]

Παρουσιάζουμε κάποιες εκδοχές της:

bGB[x_, t_] := 2x;
sigmaGB[x_, t_] := 0.1x;

EulerMaruyamaPlot[1, bGB, sigmaGB, 0.01, 3, 10,"Γεωμετρική κίνηση Brown με r=2 και c=0.1"]

bGB2[x_, t_] := 2x;
sigmaGB2[x_, t_] := 3x;

EulerMaruyamaPlot[1, bGB2, sigmaGB2, 0.01, 3, 10,"Γεωμετρική κίνηση Brown με r=2 και c=3"]
Plot
Plot

Ο αναγνώστης αξίζει να παρατηρήσει τη διαφορά των παραπάνω μοντελοποιήσεων. Σύμφωνα με το [1] (σελ. 385):

Μια κατάσταση που δείχνει πώς ο παράγοντας τυχαιότητας μπορεί να καταστρέψει ένα ισχυρά αυξητικό μοντέλο, όπως το εκθετικό.

Για τον προσδιορισμό της $X_t$ θα εφαρμόσουμε το λήμμα Itô για $G(x,t)=\ln x$. Οπότε έχουμε:

\[\begin{align*} d\ln X_t&=\left(rX_t\dfrac{1}{X_t}-\dfrac{c^2X_t^2}{2}\dfrac{1}{X_t^2}\right)dt+cX_t\dfrac{1}{X_t}dW_t\\ &=\left(r-\dfrac{c^2}{2}\right)dt+cdW_t \end{align*}\]

Ολοκληρώνοντας τις δύο μεριές από $0$ έως $t$ έχουμε:

\[ \ln X_t-\ln X_0=\left(r-\dfrac{c^2}{2}\right)t+cW_t. \]

Άρα:

\[ X_t=X_0\exp\left(\left(r-\dfrac{c^2}{2}\right)t+cW_t\right). \]

Λογιστικό μοντέλο

Εδώ θα εξετάσουμε το γνωστό μας λογιστικό πληθυσμιακό μοντέλο, αναδιαμορφωμένο ώστε να φιλοξενεί έναν παράγοντα τυχαιότητας (βλ. [1] σελ. 386-387):

\[ dX_t=rX_t\left(1-\dfrac{X_t}{K}\right)dt+cX_t dW_t \]

bLOGISTIC[x_, t_] := 5*x(1-x/6);
sigmaLOGISTIC[x_, t_] := 0.5x;

EulerMaruyamaPlot[1, bLOGISTIC, sigmaLOGISTIC, 0.01, 6, 10, "Στοχαστικό Λογιστικό Μοντέλο με r=5, K=6 και c=0.5"]
Plot

Η λύση που δίνει το [1] (σελ. 387) είναι η:

\[ X_t= \dfrac{x_0 \exp\left(\left(r-\frac{c^2}{2}\right)t+cW_t\right)}{1+x_0 \frac{r}{K}\int_0^t \exp\left(\left(r-\frac{c^2}{2}\right)\tau+cW_\tau\right)d\tau} \]

Αλλά, ας το εξαιτάσουμε καλύτερα.

⚠️ Warning
⚠️ ΥΠΟ ΚΑΤΑΣΚΕΥΗ ⚠️

Διαδικασία Ornstein-Uhlenbeck

Γνωρίσαμε τη διαδικασία Ornstein-Uhlenbeck σε προηγούμενη ενότητα, όπου βρήκαμε τη μέση τιμή και τη διασπορά της. Τώρα θα προσπαθήσουμε να βρούμε και την ίδια την $X_t$.

\[ dX_t=\theta\cdot(\mu-X_t)dt+\sigma dW_t, \ \alpha,\sigma>0. \]

Για να τη λύσουμε, θα ακολουθήσουμε διαδικασία που μοιάζει με την παραγοντική ολοκλήρωση (βλ. [4] σελ. 139-140) και εφαρμόζεται γενικότερα στις Σ.Δ.Ε. της μορφής:

\[ dX_t=(\lambda(t) X_t+\theta(t))dt+\sigma(t) dW_t \]

Αρχικά πολλαπλασιάζουμε και τις δύο μεριές της εξίσωσης με $e^{\theta t}$, οπότε έχουμε:

\[\begin{align*} & \ e^{\theta t}dX_t=e^{\theta t}\theta\cdot(\mu-X_t)dt+e^{\theta t}\sigma dW_t\\ \Leftrightarrow & \ e^{\theta t}dX_t+(e^{\theta t})'(X_t-\mu)dt=e^{\theta t}\sigma dW_t\\ \Leftrightarrow & \ e^{\theta t}d(X_t-\mu)+(e^{\theta t})'(X_t-\mu)dt=e^{\theta t}\sigma dW_t\\ \Leftrightarrow & \ d\left(e^{\theta t}(X_t-\mu)\right)=e^{\theta t}\sigma dW_t \end{align*}\]

Ολοκληρώνοντας και τις δύο μεριές από $0$ έως $t$ έχουμε:

\[ e^{\theta t}(X_t-\mu)+\mu-x_0=\int_0^t e^{\theta s}\sigma dW_s. \]

Άρα:

\[ X_t=\mu+\left(x_0-\mu+\int_0^t e^{\theta s}\sigma dW_s\right)e^{-\theta t} \]

Τεστ Feller για εκρήξεις

Η μελέτη μας θα γίνει σε ένα διάστημα:

\[ I=(\ell,r)\text{, με }-\infty\leq \ell, r\leq+\infty \]

Συγκεκριμένα, εν προκειμένω θα ασχοληθούμε με Σ.Δ.Ε. του τύπου:

\[ dX_t = b(X_t)dt+\sigma(X_t)dW_t \]

Όπου $b, \sigma : I\to \mathbb{R}$ (βλ. [2] σελ. 329 και 342-343). Ο αναγνώστης ας παρατηρήσει ότι οι συναρτήσεις $b$ και $\sigma$ δεν έχουν πλέον δύο ορίσματα, αλλά ένα.

Θα μάς απασχολήσει ο χρόνος όπου η διαδικασία μας ($X_t$) χρειάζεται για να ξεφύγει εκτός του προκαθορισμένου διαστήματος $I$. Πιο φορμαλιστικά (βλ. [2] σελ. 343):

\[ S=\inf\{t\geq 0: X_t\notin I\} \]

Αρχικά ορίζουμε (βλ. [2] σελ. 339) τη scale function:

\[ p(x)=\int_{c}^{x}\exp\left(-2\int_{c}^{\xi}\dfrac{b(\zeta)}{\sigma^2(\zeta)}d\zeta\right)d\xi \]

με $c\in\mathbb{R}$ αυθαίρετη σταθερά και υπό τις συνθήκες (βλ. [2] σελ. 343):

\[ \forall x \in I \left(\sigma^2 (x)>0\right) \]

και:

\[ \forall x \in I , \exists \varepsilon>0\text{ τέτοιο, ώστε }\int_{x-\varepsilon}^{x+\varepsilon}\dfrac{1+|b(y)|}{\sigma^2(y)}dy<+\infty \]

Ακολούθως ορίζουμε (βλ. [2] σελ. 347):

\[ v(x)=\int_{c}^{x}p'(y)\left(\int_{c}^{y}\dfrac{2dz}{p'(z)\sigma^2 (z)}\right)dy. \]

Το τεστ Feller (βλ. [2] σελ. 348) μάς λέει αν είναι πιθανό να έχουμε έκρηξη (με την έννοια ότι βγαίνει η $X_t$ εκτός του $I$), δηλαδή αν η πιθανότητα $\mathbb{P}(S<+\infty)$ είναι μη-μηδενική. Έτσι, δεν έχουμε έκρηξη (ήτοι $\mathbb{P}(S<+\infty)=0$), όταν:

\[ v(\ell+)=v(r-)=\infty, \]

αλλιώς έχουμε έκρηξη.

Αλλά, ας το δούμε και με ένα παράδειγμα.

⚠️ Warning
⚠️ ΥΠΟ ΚΑΤΑΣΚΕΥΗ ⚠️

Βιβλιογραφία

  1. Allen, Linda J. S. 2010. An Introduction to Stochastic Processes with Applications to Biology. CRC Press.
  2. Karatzas, Ioannis, and Steven E. Shreve. 1991. Brownian Motion and Stochastic Calculus. Springer.
  3. Ross, Sheldon M. 2019. Introduction to Probability Models. Academic Press.
  4. Steele, J. Michael. 2001. Stochastic Calculus and Financial Applications. Springer.
  5. Tsay, Ruey S. 2010. Ανάλυση Χρονοσειρών. Εκδόσεις Τζιόλα.
ἐπιψαύσεις - υλοποίηση με WLJS