\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{lscape}
\textwidth=7in  \textheight=15in
\hoffset=-1in \voffset=-1.25in
%\pagestyle{empty}
\newcommand{\boxspace}{\parbox{1.42in}}
\newcommand{\hhspace}{\hspace{0.02in}}
%
\newcommand{\bboxspace}{\parbox{2.0in}}
\newcommand{\hhhspace}{\hspace{0.05in}}
\input{tcilatex}
\begin{document}
\begin{flushright}
Numerical Methods of SDEs-Dr. Ra$\acute{\textrm{u}}$l Tempone\\
Homework\:$\#3$\\
\end{flushright}

\begin{description}
\item[Exercise 1]  \textit{Formulate and motivate a forward Euler method for approximation of the Stratonovich SDE}

\begin{equation}\label{eq:StatSDE}
dX_t = a(t,X_t)dt + b(t,X_t) \circ dW_t 
\end{equation}
\vspace{.25in}
\\


\textbf{First approach:}  (Incorrect approach!)  We approximated $b(t,X_t)$ by the midpoint rule:


\begin{equation}\label{eq:SDEw}
X_{t_{n+1}}=X_{n} + a(t_{n},X_{t_n})\Delta t_n + b\left(\frac {1}{2} (t_{n+1}+t_n), \frac {1}{2} (X_{t_{n+1}}+X_{t_n})\right) \Delta W_n
\end{equation}


\underline{What is wrong with approach (\ref{eq:SDEw})?}
\vspace{.25in}
\\
Consider the following process as an example to illustrate why formulation (\ref{eq:SDEw}) was incorrect.

\begin{equation}\label{eq:ex}
dX(t)= bX_tdW(t)  \hspace{.5cm} \textrm{with a constant} \, b,
\end{equation}

Using (\ref{eq:SDEw}) to approximate (\ref{eq:ex}) we obtain:

\begin{eqnarray}
X_{t_{n+1}}&= &X_{t_{n}}  + b \left(\dfrac {1}{2} (X_{t_{n+1}}+X_{t_n})\right) \Delta W_n \nonumber
\\
X_{t_{n+1}} (1-\dfrac {1}{2} b \Delta W_n) & = & X_{t_{n}} (1 + \dfrac {1}{2} b \Delta W_n) \nonumber
\\
X_{t_{n+1}} & =& X_{t_{n}}\left[\dfrac { \left(1 + \dfrac {1}{2} b \Delta W_n\right)} {\left(1-\dfrac {1}{2} b \Delta W_n\right)}\right]
\end{eqnarray}


\textbf{Consider the computation  of $E(X^2_1)$}

\begin{equation}\label{eq:int1}
E(X_1^2) = \dfrac{(X_0)^2}{\sqrt{2\pi\Delta t}} \int^{\infty}_{-\infty} \dfrac { (1 + \frac {1}{2} b z)^2} {(1-\frac {1}{2} bz )^2}\, e^{\frac{-z^2}{2 \Delta t}} dz 
\end{equation}

where we have used the fact that $E(f(x))=\int(f(x)dp(x))$ where $p(x)$ is the probability density of $x$.
\vspace{.25in}
\\
From (\ref{eq:int1}), consider a similiar integral:

\begin{equation}\label{eq:conv}
\int_{-\infty}^{\infty} \frac {(1+x)^2}{(1-x)^2}\,e^{-x^2} dx 
\end{equation}

However, the integral in (\ref{eq:conv}) does not converge.

\newpage


\textbf{Formulation:}

Rewrite (\ref{eq:StatSDE}), so that may discretize using the forward Euler method.

\begin{equation}\label{eq:StatSDE1}
X_t = X_0 + \int^t_0 a(s, X(s))ds + \dfrac{1}{2}\int^t_0 b(s,X_s) b_x(s,X_s) ds + \int^t_0 b(s,X_s)dW_s
\end{equation}
\vspace{.25cm}
\\

From (\ref{eq:StatSDE1}), considering one interval, $[t_{n},t_{n+1}]$

\begin{equation}\label{eq:StatSDE2}
X_{n+1} = X_n + \int^{t_{n+1}}_{t_{n}} a(x, X(s))ds + \dfrac{1}{2}\int^{t_{n+1}}_{t_{n}} b(s,X_s)b_x(s,X_s) ds+ \int^{t_{n+1}}_{t_{n}} b(s,X_s)dW_s
\end{equation}

Discretization of (\ref{eq:StatSDE2}) to obtain the following:

\begin{equation}\label{eq:StatSDE3}
X_{t_{n+1}} = X_{t_n} +  a(t_n, X_{t_n})\Delta t_n + \dfrac{1}{2} b(t_n, X_{t_n})\dfrac{\partial b(t_n,X_{t_n})}{\partial x}\Delta t_n + b\Delta t_n \Delta W_{t_{n}}
\end{equation}






%The formulation (\ref{eq:StatSDE1}), was chosen to approximate the Stratonovich SDE for the following reasons:


\newpage

\item[Exercise 2]  \textit{Consider the deterministic differential equation}

\begin{eqnarray}\label{eq:DDE2}
dZ_t &=& a(Z_t)dt\\
Z(0) &=& x_{0}  \hspace{.75in} 0\leq t \leq T, \nonumber
\end{eqnarray}

\textit{and consider a peturbation of} (\ref{eq:DDE2}), \textit{the It$\hat{\textrm{o}}$ stochastic differential equation}

\begin{eqnarray}\label{eq:itoDE2}
dX_t &=& a(X_t)dt + bdW_t\\
X(0) &=& x_{0}  \hspace{.75in} 0\leq t \leq T, \nonumber
\end{eqnarray}

\textit{where \textbf{a} is a smooth function and \textbf{b}\,$>0$ is a positive constant.  The aim of this exercise is to compare the solution of both equations.  Define then the difference:}

\begin{equation}\label{eq:err2}
e_{t}= X_{t} - Z_{t}\\
\end{equation}

\begin{description}
\item[a)]  \textit{Consider $\mathbf{a(x) = ax}$ (linear case) and compute $\mathbf{E(e_t)}$ and $\mathbf{Var(e_t)}$.}  
\\
\textit{Hint: Use It$\hat{\textrm{o}}$'s formula when necessary.}

From (\ref{eq:DDE2}), obtain the following solution

\begin{equation}\label{eq:solnDDE2}
Z_{t} = x_{0}e^{at}\\
\end{equation}

Applying the intergrating factor, define the following:

\begin{eqnarray}\label{eq:subDDE2}
Y_{t} &=& e^{-at}X_{t} \hspace{.75in} 0\leq t \leq T\\
Y_{0} &=& x_{0}   \nonumber
\end{eqnarray}

\textbf{It$\hat{\textrm{o}}$'s formula}. Let $f$ be a smooth function, assume that $X_{t}$
is the solution of an SDE, then $g_{t} = f(X_{t}, t)$ satisfies the SDE

\begin{equation}\label{eq:itoform}
df(X_{t}, t) = \dfrac{\partial f}{\partial t}dt +
\dfrac{\partial f}{\partial x}\,dX_{t} + 
\dfrac{1}{2}\dfrac{\partial^2 f}{\partial x^2}\,\sigma^{2}(X_{t}, t)dt,
\end{equation}

Applying (\ref{eq:itoform}) to (\ref{eq:subDDE2}), to obtain the following:


\begin{eqnarray}\label{eq:ySOLN}
dY_{t} &=& -ae^{-at}X_{t}dt + e^{-at}\underbrace{dX_{t}}_{\textrm{see} \: (\ref{eq:itoDE2})}  + \underbrace{\textrm{higher order terms}}_{= 0} \nonumber
\\
& & \nonumber
\\
&=& -ae^{-at}X_{t}dt + ae^{-at}X_tdt + be^{-at}dW_t  \nonumber
\\
dY_{t} &=&  be^{-at}dW_t 
\end{eqnarray}

integrating (\ref{eq:ySOLN}), will obtain:

\begin{eqnarray}\label{eq:solnX}
Y_{t} &=& Y_0 + \int^t_0 be^{-as}dW_{s}  \nonumber 
\\  
& & \nonumber
\\
e^{-at}X_{t} & = & X_0 + \int^t_0 be^{-as}dW_{s} \nonumber
\\
X_{t} & = & e^{at}X_0 + \int^t_0 be^{a(t-s)}dW_{s}
\end{eqnarray}

\newpage

Therefore, now that we have obtained an equation for $X_t$, see (\ref{eq:solnX}), revisit the definition for the difference, (\ref{eq:err2}).

\begin{eqnarray}\label{eq:erreq}
e_{t} &=& \underbrace{X_{t}}_{\textrm{see} \, (\ref{eq:solnX})} - \underbrace{Z_{t}}_{\textrm{see} \, (\ref{eq:solnDDE2})} \nonumber 
\\  
& & \nonumber
\\
e_{t} &=& \int^t_0 be^{a(t-s)}dW_{s}
\end{eqnarray}

\textbf{Solution:}  Now derived an equation for $e_{t}$, (\ref{eq:erreq}), compute $\mathbf{E(e_t)}$ and $\mathbf{Var(e_t)}$.
\vspace{.15cm}
\\ 
\begin{itemize}
	\item Compute $E(e_{t})$.  Will use, \textbf{Theorem 2.15 (ii)}, using the basic properities of It$\hat{\textrm{o}}$'s integrals.

\begin{displaymath}
\textrm{\textbf{Recall:}} \hspace{.25in} E\left[\int^t_0 f(s, \cdot)dW_{s}\right] = 0
\end{displaymath}

\begin{eqnarray}\label{eq:ans2a-1}
E(e_{t}) &=& E\left(\int^t_0 be^{a(t-s)}dW_{s}\right) \nonumber
\\
&=& 0  \hspace{.25in} \textbf{\textrm{by Thm 2.15}}
\end{eqnarray}
	
	\item Compute $Var(e_{t})$.  
	
\begin{eqnarray}\label{eq:ans2a-2a}
Var(e_{t}) &=& E[e^2_{t}] - [\underbrace{E[e_{t}])}_{\textrm{see} \: (\ref{eq:ans2a-1})}]^2\nonumber
\\
&=& E[e^2_{t}] \nonumber 
\\
&=& E\left[\left(\int^t_0 be^{a(t-s)}dW_{s}\right)^2\right] 
\end{eqnarray}

Will use, \textbf{It$\hat{\textrm{o}}$ isometry} in order to simplify (\ref{eq:ans2a-2a}).
\begin{displaymath}
\textrm{\textbf{Recall:}} \hspace{.25in} E\left[ \int^t_0 f(W_s, s)dW_{s}\right]^2 = 
\int^t_0 Ef^2(W_s, s)ds
\end{displaymath}

From above, (\ref{eq:ans2a-2a}) will become,
	
\begin{eqnarray}\label{eq:ans2a-2b}
Var(e_{t}) &=& E\left[\left(\int^t_0 be^{a(t-s)}dW_{s}\right)^2\right] \nonumber
\\
&=& E\left(\int^t_0 b^2e^{2a(t-s)}ds\right) \nonumber
\\
&=& \dfrac{b^2}{2a}\left(e^{2at}-1 \right) \nonumber
\end{eqnarray}
\end{itemize}

\newpage

\item[b.]  \textit{Assume now that $\mathbf{|a(x)-a(y)| \leq C_a|x-y|}$ with a positive constant $C_a$. 
Find bounds for the expectation, $\mathbf{E(|e_{t}|^2)}$  and use it to bound the variance $\mathbf{Var(e_t)}$.  \\
Discuss what happens as $b \rightarrow 0$.}

\begin{itemize}
	\item Find bounds for the expectation, $\mathbf{E(|e_{t}|^2)}$.
	
\begin{eqnarray}\label{eq:err2b}
e_t &=& X_t - Z_t \nonumber
\\
&=& \int^t_0(a(X_s) - a(Z_s))ds + \int^t_0 bdW_s 
\end{eqnarray}
\vspace{.15cm}
\\
Taking the absolute value, squaring both sides, and taking the expectation we will obtain the left hand solution.  
To obtain a upper bound, we will use various methods involving the triangle inequality, Cauchy-Schwarz inequality, Fubini's theorem, and Gr$\ddot{\textrm{o}}$nwall's lemma.

\begin{eqnarray}\label{eq:err2soln-1}
|e_t| &=& \left|\int^t_0 (a(X_s) - a(Z_s))ds + \int^t_0 b\,dW_s\right| \nonumber
\\
&\leq& \int^t_0 C_a|X_s - Z_s|ds + b\,|W_t| \hspace{.25in} (b>0, \textrm{by assumption})  \nonumber
\\
& = & \int^t_0 C_a|e_s|ds + b\,|W_t|  \nonumber
\\
|e_t| &\leq& \int^t_0 C_a|e_s|ds + b\,|W_t|
\end{eqnarray}	

From (\ref{eq:err2soln-1}), we will square both sides, and take the expectation to obtain:

\begin{eqnarray}\label{eq:err2soln-2}
E|e_t|^2 &\leq & E\left[\left(\int^t_0C_a|e_s|ds + b\,|W_t|\right)^2\right] \nonumber
\\
&\leq&  2E\left[\left(\int^t_0C_a|e_s|ds\right)^2\right] + 2E\left[\left(b\,|W_t|\right)^2\right]
\hspace{.25in} \textbf{\textrm{Recall:}} \:\: (a+b)^2\leq 2(a^2+b^2)   \nonumber
\\
& = & 2E\left[\left(\int^t_0 C_a|e_s|ds\right)^2\right] + 2E\left[b^2\,W^2_t\right]  \nonumber
\\
& \leq & 2E\left[\left(\int^t_0 C^2_ads\right)\left(\int^t_0 |e_s|^2ds\right)\right] + 2b^2t  \hspace{.25in} \textbf{\textrm{Recall:} Cauchy-Schwarz inequality} \nonumber
\\
& \leq & 2C^2_atE\left[\int^t_0 |e_s|^2ds\right] + 2b^2t \nonumber
\\
E|e_t|^2 & \leq & 2C^2_at\left[\int^t_0 E|e_s|^2ds\right] + 2b^2t \hspace{.25in} \textbf{\textrm{Recall:} Fubini's theorem } \:\: %\href{www.}{text}
\end{eqnarray}
%postive and continous w/r.t s

To simplify the integral in (\ref{eq:err2soln-2}), we will use Gr$\ddot{\textrm{o}}$mwall's lemma.
\vspace{.25in}
\\
\textbf{Lemma.}  Assume that there exist positive constants $A$ and $K$ such that the function $f : \mathbb{R} \rightarrow \mathbb{R}$ satisfies:
\begin{eqnarray}
f(t) &\leq& K\int^t_0 f(s)ds + A \nonumber
\\
\textrm{then} & & \nonumber
\\
 f(t) &\leq&  Ae^{Kt}
\end{eqnarray}

\newpage

Also, to use Gr$\ddot{\textrm{o}}$nwall's lemma, define the following:
\begin{equation}\label{eq:func}
f(t) = E|e_t|^2
\end{equation}


From (\ref{eq:func}), and substituting into (\ref{eq:err2soln-2}), we will obtain:

\begin{eqnarray}\label{eq:err2soln-3}
f(t) & \leq & 2C^2_at\int^t_0 f(s)ds + 2b^2t, \nonumber
\\
& \leq & 2C^2_aT\int^t_0 f(s)ds + 2b^2T
\end{eqnarray}

Let $K =2C^2_aT$ and $A = 2b^2T$, \:\: (note: $A >0$)

Then by Gr$\ddot{\textrm{o}}$nwall's lemma, obtain a bound on $\mathbf{E(|e_{t}|^2)}$.

\begin{equation}\label{eq:ANS2b}
f(t) \leq 2b^2T\,e^{2C_a^2Tt}
\end{equation}
	
\newpage
\item Using answer from part 2a), (\ref{eq:ANS2b}), obtain a bound on the variance $\mathbf{Var(e_t)}$.

\begin{eqnarray}\label{eq:var1}
\textrm{Var}(e_t) &=& E(e^2_t) - \underbrace{[E(e_t)}_{\textrm{see}\,(\ref{eq:ans2a-1})}]^2 \nonumber
\\
& \leq & E(e^2_t) \hspace{.75cm} (Note: [E(e_t)]^2 \geq 0) \nonumber
\\
\textrm{Var}(e_t) &\leq & 2b^2T\,e^{2C_a^2Tt}
\end{eqnarray}
	
\newpage
	
\item Discuss what happens as $b \rightarrow 0$.
\vspace{.25cm}
\\
From (\ref{eq:var1}), as $b\rightarrow 0$, will lead to:
\begin{equation}
0 \leq \textrm{Var}(e_t) \leq 0  \Longrightarrow \textrm{Var}(e_t) = 0
\end{equation}
	
Therefore, the stochastic term disappears as $b\rightarrow 0$ in the SDE of $X_t$.  Thus, matching the deterministic
function, $Z_t$
\end{itemize}

\newpage

\item[c)] \textit{Implement a uniform time step forward Euler discretization of the above equations taking
\\
$\mathbf{a(x) = cos(x), b=0.1,}$ and $\mathbf{T = 6}$.  
\\
Plot the sample estimator for $\mathbf{Var(e_t)}$ vs. time.
\\
Compare it with the bound obtained in part (b).  
\\
Use $\mathbf{M = 10^3}$ sample paths and different number of time steps: $\mathbf{N = 10,20, 40}$.}
\end{description}

\begin{itemize}
	\item Each graph below was graphed using Matlab, please see appendix for code.
	\item  For each figure, the deterministic function, trajectories, and the bounds were plotted.
	\vspace{.15cm}
	\\
\end{itemize}

\begin{itemize}
	\item  \textbf{Figure one:}  Number of time steps: $\mathbf{N = 10}$.  
\end{itemize}
\end{description}
	\begin{flushleft}
\pdfimage width 20cm {n10.png}
	\end{flushleft}


\newpage
\begin{itemize}
	\item  \textbf{Figure two:}  Number of time steps: $\mathbf{N = 20}$.  
\end{itemize}
	\begin{flushleft}
\pdfimage width 20cm {n20-2.png}
	\end{flushleft}


\newpage
\begin{itemize}
	\item  \textbf{Figure three:}  Number of time steps: $\mathbf{N = 40}$.  
\end{itemize}
	\begin{flushleft}
\pdfimage width 20cm {n40.png}
	\end{flushleft}


\newpage
\begin{center}
\LARGE{APPENDIX:}  Problem 2(c) Mathlab code
\end{center}
\begin{verbatim}
clf randn(round(clock),1); T = 6; N = 40; Delta_t = T/N; b=0.1; X0=0;
C1=2*(b^2); C2=8; v_error_mean=zeros(N,1); v_error_square=zeros(N,1);
v_variance=zeros(N,1); v_xdet=zeros(N,1); v_xdet2=zeros(N,1); M = 1e+3;
for s = 1:M
    X=X0;
    Z=X0;
    for j = 1:N
        Winc1 = sqrt(Delta_t)*randn;
        X=X+cos(X)*Delta_t + b*Winc1;
        Z=Z+cos(Z)*Delta_t;
        v_xdet(j,1)=v_xdet(j,1)+X;
        v_xdet2(j,1)=v_xdet2(j,1)+X^2;
        e=abs(X-Z);
        v_error_mean(j,1)=v_error_mean(j,1)+e;
        e2=(X-Z)^2;
        v_error_square(j,1)=v_error_square(j,1)+e2;
    end end X=X0; X1=X0; Z=X0; for j=1:N
        v_xdet(j,1)=v_xdet(j,1)/M;
        v_xdet2(j,1)=v_xdet2(j,1)/M;
        v_variance_det(j,1)=v_xdet2(j,1)-((v_xdet(j,1))^2);
        Winc1 = sqrt(Delta_t)*randn;
        X=X+cos(X)*Delta_t + b*Winc1;
        v_notmean(j,1) = X;
        Winc2 = sqrt(Delta_t)*randn;
        X1=X1+cos(X1)*Delta_t + b*Winc2;
        v_notmean2(j,1) = X1;
        v_error_mean(j,1)=v_error_mean(j,1)/M;
        v_error_square(j,1)=v_error_square(j,1)/M;
        v_variance(j,1)=v_error_square(j,1)-((v_error_mean(j,1))^2);
        Z=Z+cos(Z)*Delta_t;
        v_deterministic(j,1)=Z;
        v_time=linspace(0,T,N);      
        v_minus_2d(j,1)=v_xdet(j,1)-2*sqrt(v_variance_det(j,1));
        v_plus_2d(j,1)=v_xdet(j,1)+2*sqrt(v_variance_det(j,1)); end
title('N = 10'); subplot(2,1,1); plot(v_time, v_variance,'r');
legend('Variance(e(t))'); xlabel('Time (in seconds)'); subplot(2,1,2);
plot(v_time, v_deterministic, 'k', v_time, v_xdet,'b', v_time, v_notmean,
'm+--', v_time, v_notmean2,'gx--'); hold on; plot(v_time, v_minus_2d,'r',
v_time, v_plus_2d,'r','LineWidth',3); legend('Z','Avg(X)','Trajectory
1','Trajectory 2','Avg(X) +/- 2 SD');
xlabel('Time (in seconds)');
\end{verbatim}




%end letters




%ends ex #
\end{document}