\documentclass[notitlepage]{article} \usepackage{../../Math556} \usepackage{listings} \usepackage{numprint} \usepackage{enumerate} \usepackage{multirow} \usepackage{bbm} \usepackage{bbold} \usepackage{amsfonts,amsmath,amssymb} \setlength{\parindent}{0pt} \def\E{\Expect} \lstloadlanguages{R} \definecolor{keywordcolor}{rgb}{0,0.6,0.6} \definecolor{delimcolor}{rgb}{0.461,0.039,0.102} \definecolor{Rcommentcolor}{rgb}{0.101,0.043,0.432} \lstdefinestyle{Rsettings}{ basicstyle=\ttfamily, breaklines=true, showstringspaces=false, keywords={if, else, function, theFunction, tmp}, % Write as many keywords otherkeywords={}, commentstyle=\itshape\color{Rcommentcolor}, keywordstyle=\color{keywordcolor}, moredelim=[s][\color{delimcolor}]{"}{"}, } \lstset{basicstyle=\ttfamily, numbers=none, literate={~} {$\sim$}{2}} \def\htitle{\textsclarge{Some examples of characteristic Functions}} \begin{document} <>= library(knitr) # global chunk options opts_chunk$set(cache=TRUE, autodep=TRUE) options(scipen=999) options(repos=c(CRAN="https://cloud.r-project.org/")) inline_hook <- function (x) { if (is.numeric(x)) { # ifelse does a vectorized comparison # If integer, print without decimal; otherwise print two places res <- ifelse(x == round(x), sprintf("%.6f", x), sprintf("%.6f", x) ) paste(res, collapse = ", ") } } knit_hooks$set(inline = inline_hook) @ \begin{center}\coursetitle\end{center} \begin{center}\htitle\end{center} The characteristic function for a random variable $X$ with pmf/pdf $f_X$ is defined for $t \in \R$ as \begin{align*} \varphi_X(t) & = \Expect_X[e^{itX}] = \Expect_X[\cos(tX) + i \sin(tX)] = \Expect_X[\cos(tX)] + i \Expect_X[\sin(tX)] \\[6pt] & = \int_{-\infty}^\infty e^{itx} \ dF_X(x) = \int_{-\infty}^\infty \cos (tx) \ dF_X(x) + i \int_{-\infty}^\infty \sin (tx) \ dF_X(x) \end{align*} using the $dF_X(x)$ notation, where as usual the `integral' is a sum in the discrete case. For an arbitrary discrete distribution, the pmf takes the form \[ f_X(x) = \sum_{j=1}^\infty \omega_j \One_{\{x_j\}}(x) \] where $x_1,x_2,\ldots$ are a countable collection of (distinct) real values, and $0 \leq \omega_j \leq 1$ and $\sum\limits_{j=1}^\infty \omega_j = 1$. We have by direct calculation that \[ \varphi_X(t) = \sum_{j=1}^\infty \omega_j \exp\{i t x_j\} = \sum_{j=1}^\infty \omega_j (\cos(t x_j) + i \sin(t x_j)) = \sum_{j=1}^\infty \omega_j \cos(t x_j) + i \sum_{j=1}^\infty \omega_j \sin(t x_j). \] Thus \begin{align*} |\varphi_X(t)|^2 & = \left(\sum_{j=1}^\infty \omega_j \cos(t x_j) \right)^2 + \left(\sum_{j=1}^\infty \omega_j \sin(t x_j) \right)^2 \\[6pt] & = \sum_{j=1}^\infty \omega_j^2 \cos^2(t x_j) + 2 \sum_{j=1}^\infty \sum_{k=j+1}^\infty \omega_j \omega_k \cos(t x_j) \cos(t x_k) \\[6pt] & \qquad + \sum_{j=1}^\infty \omega_j^2 \sin^2(t x_j) + 2 \sum_{j=1}^\infty \sum_{k=j+1}^\infty \omega_j \omega_k \sin(t x_j) \sin(t x_k) \\[6pt] & = \sum_{j=1}^\infty \omega_j^2 + 2 \sum_{j=1}^\infty \sum_{k=j+1}^\infty \omega_j \omega_k \cos(t (x_j - x_k)) \end{align*} It can be shown that in this case that $\limsup\limits_{|t| \longrightarrow \infty} |\varphi_X(t)| = 1$; in fact, as $|\varphi_X(t)| \leq 1$, it suffices to show that \[ \limsup_{|t| \longrightarrow \infty} \text{Re} \ \varphi_X(t) \geq 1. \] where Re indicates the real part of the complex-valued function. Here \[ \text{Re} \ \varphi_X(t) = \sum_{j=1}^\infty \omega_j \cos(t x_j). \] \bigskip \textbf{Example:} $X \sim Poisson(\lambda)$ \[ \varphi_X(t) = \exp\{\lambda(e^{it}-1)\} \] This function has a real and imaginary part. We have that \[ \varphi_X(t) = \exp\{\lambda (\cos(t)-1) + i \lambda \sin(t)\} = \exp\{\lambda (\cos(t)-1)\}(\cos(\lambda \sin(t))+ i \sin(\lambda\sin(t))) \] so that \[ |\varphi_X(t)| = \exp\{\lambda (\cos(t)-1)\}. \] Therefore the function is periodic with period $2 \pi$, and \[ \limsup_{|t| \longrightarrow \infty} |\varphi_X(t)| = 1, \] as is typical for a discrete distribution. \pagebreak For $\lambda=2$, we have the following plots <>= tv<-seq(-10,10,by=0.01) lambda<-2 varphi<-exp(lambda*(exp(1i*tv)-1)) par(mar=c(4,4,2,0)) plot(tv,Re(varphi),type='l',axes=FALSE, ylim=range(-1,2),xlab='t',ylab=' ') #Real part axis(side=1,at=seq(-10,10,5)) axis(side=2,at=c(-1,0,1)) box() lines(tv,Im(varphi),col='red') #Imaginary part lines(tv,Mod(varphi),col='blue') #Modulus legend(-10,2,c(expression(Re(varphi[X](t))), expression(Im(varphi[X](t))), expression(Mod(varphi[X](t)))),col=c('black','red','blue'),lty=1) title(expression(varphi[X](t))) @ \pagebreak For $\lambda=10$: <>= lambda<-10 varphi<-exp(lambda*(exp(1i*tv)-1)) par(mar=c(4,4,2,0)) plot(tv,Re(varphi),type='l',axes=FALSE, ylim=range(-1,2),xlab='t',ylab=' ') #Real part axis(side=1,at=seq(-10,10,5)) axis(side=2,at=c(-1,0,1)) box() lines(tv,Im(varphi),col='red') #Imaginary part lines(tv,Mod(varphi),col='blue') #Modulus legend(-10,2,c(expression(Re(varphi[X](t))), expression(Im(varphi[X](t))), expression(Mod(varphi[X](t)))),col=c('black','red','blue'),lty=1) title(expression(varphi[X](t))) @ \pagebreak \textbf{Example:} $X \sim Exponential(\lambda)$ \[ \varphi_X(t) = \frac{\lambda}{\lambda - it} = \frac{1}{1-i t \mu} = \frac{1}{1+\mu^2 t^2} + i \frac{\mu t}{1+\mu^2 t^2} \] where $\mu = 1/\lambda$. Here \[ |\varphi_X(t)| = \frac{1}{\sqrt{1+ \mu^2 t^2}} \longrightarrow 0 \quad \text{as } |t| \longrightarrow \infty. \] For $\lambda=1$ <>= lambda<-1;mu<-1/lambda varphi<-1/(1-1i*tv*mu) par(mar=c(4,4,2,0)) plot(tv,Re(varphi),type='l',axes=FALSE, ylim=range(-1,2),xlab='t',ylab=' ') #Real part axis(side=1,at=seq(-10,10,5));axis(side=2,at=c(-1,0,1));box() lines(tv,Im(varphi),col='red') #Imaginary part lines(tv,Mod(varphi),col='blue') #Modulus legend(-10,2,c(expression(Re(varphi[X](t))), expression(Im(varphi[X](t))), expression(Mod(varphi[X](t)))),col=c('black','red','blue'),lty=1) title(expression(varphi[X](t))) @ For $\lambda=5$ <>= lambda<-5;mu<-1/lambda varphi<-1/(1-1i*tv*mu) par(mar=c(4,4,2,0)) plot(tv,Re(varphi),type='l',axes=FALSE, ylim=range(-1,2),xlab='t',ylab=' ') #Real part axis(side=1,at=seq(-10,10,5)) axis(side=2,at=c(-1,0,1));box() lines(tv,Im(varphi),col='red') #Imaginary part lines(tv,Mod(varphi),col='blue') #Modulus legend(-10,2,c(expression(Re(varphi[X](t))), expression(Im(varphi[X](t))), expression(Mod(varphi[X](t)))),col=c('black','red','blue'),lty=1) title(expression(varphi[X](t))) @ \pagebreak \textbf{Example:} Mixed distribution: suppose \[ F_X(x) = \frac{1}{4} \One_{[0,\infty)}(x) + \frac{3}{4} \One_{(0,\infty)}(x) ( 1 - e^{-x}) \] which places a mass of probability of size 0.25 at $x=0$, and then adds an $Exponential(1)$ distribution with probability 0.75 to contribute probability for $x > 0$. <>= x<-seq(0,5,by=0.01) Fx<-0.25*(x>=0)+0.75*(x > 0)*(1-exp(-x)) par(mar=c(4,4,2,0)) plot(x,Fx,type='l',xlim=range(-1,5),ylim=range(0,1),ylab=expression(F[X](x))) lines(c(-1,0),c(0,0)) points(0,0.25,cex=1.0,pch=19) @ Then \[ \varphi_X(t) = \int_{-\infty}^\infty e^{itx} dF_X(x) = \frac{1}{4} e^{it0} + \frac{3}{4} \int_0^\infty e^{itx} e^{-x} \ dx = \frac{1}{4} + \frac{3}{4} \frac{1}{1-it} \] so therefore \[ \varphi_X(t) = \frac{1}{4} + \frac{3}{4} \frac{1}{1+t^2} + i \frac{3}{4} \frac{t}{1+t^2} = \frac{4+t^2}{4(1+t^2)} + i \frac{3t}{4(1+t^2)} \] Here \[ |\varphi_X(t)|^2 = \frac{(4+t^2)^2}{16(1+t^2)^2} + \frac{9t^2}{16 (1+t^2)^2} = \frac{t^4 + 17t^2 + 16}{16 (1+t^2)^2} = \frac{t^2 + 16}{16(1+t^2)} \] so hence \[ \lim_{|t| \longrightarrow \infty} |\varphi_X(t)| = \frac{1}{4}. \] Thus the limiting behaviour of this mixed-type distribution is not the same as for the purely discrete or purely (absolutely) continuous type. \pagebreak <>= varphi<-0.25+0.75/(1-1i*tv) par(mar=c(4,4,2,0)) plot(tv,Re(varphi),type='l',axes=FALSE, ylim=range(-1,2),xlab='t',ylab=' ') #Real part axis(side=1,at=seq(-10,10,5));axis(side=2,at=c(-1,0,1));box() lines(tv,Im(varphi),col='red') #Imaginary part lines(tv,Mod(varphi),col='blue') #Modulus legend(-10,2,c(expression(Re(varphi[X](t))), expression(Im(varphi[X](t))), expression(Mod(varphi[X](t))), 'Limiting value 1/4'),col=c('black','red','blue'),lty=c(1,1,1,2)) title(expression(varphi[X](t))) abline(h=0.25,lty=2) @ \pagebreak \textbf{Numerical Inversion:} We may also attempt numerical inversion of a given cf to get at the pmf or pdf. Suppose \[ \varphi_X(t) = \exp\{-|t|\} \] which if the cf for the standard Cauchy distribution. Now the relevant inversion formula is \[ f_X(x) = \frac{1}{2 \pi} \int_{-\infty}^\infty e^{-itx} \varphi_X(t) \ dt = \frac{1}{2 \pi} \int_{-\infty}^\infty e^{-itx} \exp\{-|t|\} \ dt \] which we may approximate by \[ \widehat f (x) = \frac{1}{2 \pi h} \sum_{j=1}^N e^{-it_j x} \exp\{-|t_j|\} \] for a suitable grid $t_1,\ldots,t_N$ with $t_{j+1}-t_j = h$, and $N$ large enough. <>= par(mar=c(4,4,2,0)) x<-seq(-10,10,by=0.01) fx.hat<-x*0 for(i in 1:length(x)){ fx.hat[i]<-fx.hat[i]+Re(sum(exp(-1i*tv*x[i])*exp(-abs(tv)))) } fx.hat<-0.5*fx.hat*(tv[2]-tv[1])/pi plot(x,fx.hat,type='l',xlab='x',ylab=expression(hat(f)(x))) lines(x,dt(x,1),col='blue') @ If \[ \varphi_X(t) = \exp\{-|t|^\alpha\} \] for $0 < \alpha \leq 2$, the corresponding distribution is still absolutely continuous, but the form of the pdf is not standard. We can use the numerical method. <>= par(mar=c(4,4,2,0)) x<-seq(-10,10,by=0.01) tv<-seq(-50,50,by=0.01) alpha<-1.75 fx.hat<-x*0 for(i in 1:length(x)){ fx.hat[i]<-fx.hat[i]+Re(sum(exp(-1i*tv*x[i])*exp(-abs(tv)^alpha))) } fx.hat<-0.5*fx.hat*(tv[2]-tv[1])/pi plot(x,fx.hat,type='l',xlab='x',ylab=expression(hat(f)(x))) @ \end{document}