APROXIMACION DE PROCESOS DE DIFUSION148.206.53.84/tesiuami/reportesok/UAMR0823.pdf · Entonces, la...

Vol. I 1

APROXIMACION DE PROCESOS

DE DIFUSION

Colectivo de Autores*

No. 2

*Bazúa Rueda , Luis Fe1 ipe

3

133 3

Hernández Castaños , Di ego Bri ci o

Torres Cházaro, Adolfo

Velasco Hernández, Jorge X. Villa Arce, Enrique

2,5

194 Este trabajo fue parcialmente apoyado por DGICSA.SEP. a través del Convenio

C88-01-0084 que tiene esta Institución

con el área de Análisis Aplicado de la

UAM- I

( 1 ) Maestría en Estadística e Investigación de Operaciones. Colegio de Ciencias y Humanidades-IIMAS, UNAM.

(2) Maestría en Matemáticas, UAM-Iztapalapa.

(3) Departamento de Matemáticas, UAM-Iztapalapa.

(4) Escuela de Ciencias Fisicomatemáticas Universidad Autónoma de Sinaloa.

(5) Departamento El Hqmbre y su Ambiente, UAM-Xochimilco.

INTRODUCCION

I . PROCESOS ESTOCASTICOS Y ECUACIONES DIFERENCIALES 1.1 So luc ión P robab i l í s t i ca de Ecuaciones Diferenciales

1.2 E l movimiento browniano

1.3 Difusiones en una dimensión 1.4 Soluciones de problemas de Ecuaciones Diferenciales Parciales

I I . SOLUCION NUMERICA DE ECUACIONES DIFERENCIALES ESTOCASTICAS 2.1 Introducción 2.2 Efecto de l a d i s c re t i zac i ón del tiempo 2.3 Efecto de l a aproximación del MB 2,4 Resultados computacionales

I I I . GENERACION 3E i'lOVIMIENT0 BROWNIAN0

3.1 Introducción 3.2 Aproximación mediante caminatas aleatori as 3.3 Simulación de !dB por caminatas aleatorias

3.4 Aproximaciones mediante PTU.

3.5 Resul tados computacional es

I V . APROXIMACION MEDIANTE CMTC

4.1 Introducción 4.2 Método de D iscret izac ión

4,3 Convergencia del Método

4.4 Resultados Computacionales

V . INTERPRETACION BIOLOGICA DE LOS RESULTADOS

APENDICE A Programas para generar MB APENDICE B Programas para generar difusiones

BIBLIOGRAFIA

I NTRODUCC I ON

En es te t rabajo presentamos algunos métodos para la construcción de aproximaciones a procesos de difusión, ilustrándolos con resultados computacionales. Aun cuando contar con ta les técnicas computacionales resul ta de uti 1 idad en una gran variedad de campos , hemos motivado l a - discusión en el capítulo 1 en términos de l a necesidad de evaluar soluciones de problemas de valores a la f ron tera (PVF) y problemas de valores in ic ia les y a la f rontera (PVIF) asociados a operadores diferenciales de segundo orden. A grandes rasgos, l a organización del trabajo es 1 a sigui ente:

Capítulo 1. Se comienza por dar algunos ejemplos que i lus t ran l a representación estocásti ca de soluciones de los problemas antes menciona - dos, todo en forma casuística y con hnimo exploratorio. S in excepción, los problemas abordados se refieren al operados Laplaciano, en dimensio- nes 1 y 2. Enseguida , se da una breve introducción informal al movimien - to browniano, tanto desde e l punto de v i s t a f í s i c o como en referencia al modelo de Wiener. Finalmente, se presentan los procesos de difusión en forma t a l que su conexión con los operadores diferenciales de segundo - orden resulta evidente, siguiendo el enfoque del "problema de la martina gala" de Stroock y Varadhan (1979). En e s t a última sección se l lega a - 1 a representación de los procesos de d i fusión como soluciones de ecuacio - nes diferenciales estocásticas ( E D E ) , haciendo evidente su carácter marko - ffiano, del cual se deducen las representaciones estocásticas antes men - ci onadas .

-

Capítulo 2. Se parte de la representación de las difusiones como soluciones de EDE's con e l f i n de diseñar algoritmos computacionales que produzcan buenas aproximaciones de las t rayector ias de dichos procesos. Con tal propósito se retoman algunos resultados de Rllmel i n (1982), traba - j o en e l que se investiga el efecto de la discret ización del tiempo en - 1 a EDE medi ante la fami 1 i a de métodos de Runge Kutta, haciendo ver l a - necesidad de introducir un término de corrección apropiado para lograr -

una buena convergencia. Dicho trabajo supone l a d i s pon i b i l i d ad de tra-

yectorias del proceso de Wiener, l o cual en l a p rác t i ca r e su l t a menos que rea l i s t a , ya que só lo pueden generarse aproximaciones más o menos suaves

de dichas trayectorias El efecto de ta l subst i tuc ión se toma 1 uego en

cuenta, s iguiendo los trabajos de Janssen (1984). Finalmente, se presen - t an l o s r e su l t ado s de la implantación computacional de las técnicas ante - r i o re s .

Capítulo 3. Queda establecida en el capítulo 2 l a necesidad de

contar con buenas aproximaciones de t rayector ias del proceso de Wiener

(MB) para generar las de procesos de d i fu s ión . En este capítulo se obtiene un paquete computacional para la generación de trayector ias del - MB, basándose en dos ideas teór icas ya establec idas: la pos ib i l idad de

aproximar el MB mediante a) caminatas aleatorias LKnight (1981)], b) pro - ceses de transporte uniforme (PTU) [Griego e t a l (1971) I , [Gorostiza y

Griego (1980) l . Por ambos caminos se obtiene convergencia uniforme en

interva los compactos, con probabilidad 1, l o cual es más que su f i c iente para garant izar la bondad de l a aproximación en lo s cá lcu los que se moti - van en el capítulo 1. Se incluyen resultados de l a simulación, que i n - cluyen una comparación de ambos métodos para generar aproximaciones del

MB.

Capítulo 4. En este ú l t imo capí tu lo se examina una ruta alternativa para aproximar trayectorias de d i fus iones, misma que genera aproximaciones discretas que convergen débilmente a la d ifus ión; esta conver-

genc ia es jus to la que se necesita para garanti zar la bondad de l a s a p r g ximaciones a l a s so luc iones e s tocás t i cas de l o s problemas que s i r ven de motivación en el capítulo 1. El método se basa en la d i sc ret i zac ión de l

espacio mediante una técnica dada en [Di Masi/Runggaldier (1980)1, l a cual da lugar a una cadena de Markoff a tiempo continuo (CMTC). A su - vez, esta técnica s igue.de las consideradas por [Kushner (1977)], basa - das en l a s ub s t i t uc i ón de l a d i f u s i ó n p o r una cadena de Markoff a tiempo discreto. En el capítulo 10 de [Hernández (1985)l hemos considerado las aplicaciones de este ti po de aproximaciones discretas en el d iseño de -

http://sigue.de

f i l t r o s para la estimación recursiva del estado de sistemas dinámicos esto cást icos , lo que contribuye a l a motivación para considerar este tipo de - técnicas. Se incluyen resultados de l a implantación computacional de estas técni cas.

Para terminar y a 1 a vez redondear esta Introducción , i ndicaremos l a manera en que las trayectorias generadas en los capí tulos 2, 3 y 4 pue - den u t i l i za r se para aproximar valores de la solución de W F ' s y PVIF's, - mediante las respectivas técnicas de Monte Carlo. Para e l l o y, sin s a l i r - nos del ámbito unidimensional , sea el operador elíptico,

1 d 2 d A:= a ( x ) - + b ( x ) - dx dx

y sea T el instante en que 1 a difusión x. asociada a A abandona e l i n - tervalo D:=(a,b) por vez primera. Finalmente, sea Px la probabilidad

según l a cual P (x =x) = 1 y denotemos Ex l a operación de tomar espe-

ranza con respecto a Px. x 0

En la sección 1.4 s e demuestra que l a fórmula

da la solución del W F

Au=O en D; u l a D = f

Asimismo, en dicha sección también se demuestra que 1 a solución del PVIF

a = Au en ( 0 , ~ ) xD a t

u (0 ,x ) = f ( x ) , x E D

Dicho PVIF resu l ta , por ejemplo, en el siguiente contexto: se tiene una barra metálica uniforme de l o n g i t u d 1 y difusividad térmica a , aislada térmicamente por todas partes. Además, se conoce el perfi 1 i n i - c ia l de temperatura a l o largo de l a barra., Se desea determinar cómo va - r í a en el tiempo dicho perfi l de temperatura. Para ello, sea u(t ,x):=tem - peratura en el p u n t o x de l a barra, t segundos después de iniciado el experimento y supongamos que u t iene e l grado requerido de suavidad, - Entonces, la ley de la conservación de la energía acoplada con l a ley de Fourier para l a conducción del calor se traduce precisamente en e s t e PVIF. Para l o s deta l les , véase [Bird e t , a l (1960)1,

De nuevo por e l método de separación de variable, obtenemos la solu - ción en l a forma de 1 a s e r i e de Fourier

W

u ( t , x ) = A. + C \ exp (-a - k nx k = l

donde 1 A. = 1s f ( s ) d s

l o 1

O Ak = ;{ f (S) cos 1 k m ds, k = l , 2 . . .

Aquí también podemos substituir las expresiones que acabamos de da r para los coeficientes de 1 a s e r i e de Fourier, i ntercambi ando 1 uego e l orden de la suma y l a i n t eg ra l , para obtener así la solución en l a forma

u ( t , x ) = j o P t ( x , s ) f ( s ) d s .

don de

f1 Pt(x,s) > O , O < S < 1 ; Pt(x,s) f(s) ds = 1.

- - J o

Razonando igual que en los ejemplos anteriores, construimos una va - r iable a leator ia Xt d is t r ibuida según

/b P(a < Xt < b ) = j Pt(x,s) ds.

a

con lo cual podemos dar la sol uci ón en 1 a forma

Se l lega as í a una expresión probabilísti ca para la solución del PVIF planteado: el perfil de temperatura en cada instante se determina promedi ando adecuadamente los valores del perfi 1 i n i cia1 o

Hay semejanzas y diferencias entre las fórmulas (1) y ( 2 ) por un lado y l a fórmula (3) . En las siguientes secciones exploraremos estas relaciones en términos del concepto de proceso estocástico. Por lo pron - t o podemos anotar ya el hecho de que pueden darse expresiones probabi l i s - t i ca s para las soluciones de PVF's y PVIF's asociados a algunos operadores diferenciales de segundo orden. Aunque aquí hemos i lustrado l o anter ior sólo para tres operadores específicos, en la sección 4 haremos ver cómo estos resultados se desprenden como coro1 ar ios de otros, mucho más generales. Por lo pronto, daremos una breve introducción al campo - de los procesos estocásticos en las s iguientes dos secciones, con objeto de tener a l a mano el material requerido para estos fines,

Agradecemos la eficiente labor mecanográfica de este t rabajo real i - zado por Beatriz Arce y Sonia Cabeza.

L A P I T U L O I

1. Soluciones probabi l í st icas de ecuaciones diferenciales

Los ejemplos que damos a cont inuac ión i lust ran las muchas interac - cienes que existen entre campos de l a Matemática aparentemente disconexos como son la Probabi 1 idad y l a s Ecuaciones en Derivadas Parciales. Más - adelante, en l a última sección de este capítulo, haremos ver cómo pueden deducirse de manera s istemática muchos otros ejemplos análogos. Todos

estos resultados pertenecen a una rama de l Aná l i s i s de gran r iqueza, cu-

yos orígenes pueden t razar se a [Kac (1951) , (1966)] , [Kakutan (1944)l y,

aún más a t rás en el tiempo, a [Courant et a l (1928) l .

EJEMPLO 1. Consideremos el problema de valores a l a f r o n t e r a (P'JF)

u" = O en (a,b) 9 u(a) = A,u(b) = B

cuya solución es

u(x) = A(&) + B ( E ) .

según se comprueba con faci l idad. Observamos que l o s dos términos entre

paréntesis son no negativos para a " < x < b, además de que suman 1; deno

temos al primero mediante el símbolo px.

Podemos pensar en el siguiente experimento: fija x, construimos

una moneda que, a l lanzar la, ar roje águ i la con probabilidad px, sol con

probabil idad 1-px: s i sale águi la, e legimos el extremo a del interva -

l o ta ,b3, e l o t ro extremo en caso contrario. Procediendo de manera más formal, sea Z una VA t a l que

P(Z = a) = px, P(Z = b) = 1-px.

y sea f l a función dada por al -> A, b l --> B. Entonces, f(Z) es una var iab le a leator ia d i screta, que toma l o s va l o re s A y B con pro-

2.

babil idad px y 1-p , respectivamente: así pues, su esperanza no es otra

que Apx + B ( l - p ), es decir , u ( x ) . En otras palabras, X

u ( x ) = Exf(Z)

Hemos incluido el subindice x para hacer exp l í c i to e l hecho de que el punto x permanece f i j o en toda 13 d i scusión,

EJEMPLO 2, Sea ahora e l W F

Au = O en D , U l a D = f

donde D es u n disco abierto de radio a en e l plano, f es - continua en el contorno de D y A denota el operador Lap1 aci ano, Basta introducir coordenadas polares (rye) y apl icar e l método de separación de variables [Courante-Hi l b e r t , vol. I(1966)j seguido de una manipulación - formal intercambiando se r i e s e integrales, para 1 legar a que

expresión que da el valor de la solución en (r ,e) como un promedio adecua - do de los valores que toma en la frontera: en efecto, resul ta que.

1 a2 - r2 . r IT a2 + r2 - 2 arcost

P ( t ) : = 2 O < t < Z I T ”

(el llamado núcleo de Poisson), de donde se ve que

además de que

3.

Esta última relación es Máximo [Levi (1980) , p. 421 o

Igual que en el ejemplo

consecuencia directa del Principio del -

anter ior , podemos construi r una V A Z que " e l i j a u n punto en a D l ' según l a l ey de probabilidad

P(a < Z < b ) =

con l o cual de nueva cueyta obtenemos que

A q u í hemos incl uído el subindice (r,e ) para indicar que l a l ey - probabil ística correspondiente cambia con el punto elegido.

Podemos observar un buen número de semejanzas entre los dos ejemplos anteriores. Por ejemplo, en ambos la expresión probabilist ica de - la solución es del tipo "esperanza de los valores de l a condición a l a frontera calculada en puntos elegidos aleatoriamente sobre l a misma", - con una ley probabil íst i ca dis t inta para cada punto del dominio en cues- tión. Por otro lado, hay que reconocer que ambos ejemplos tratan de sen- dos PVF's asociados a operadores diferenciales de segundo orden: es to puede muy bien estar en la ra íz de l a s semejanzas observadas entre (1) y ( 2 ) . Veamos a continuación un ejemplo de naturaleza dis t inta , en el que se resuelve un problema de valores iniciales y a la f rontera (PVIF) asociado a un operador 1 ineal paraból i co.

EJEMPLO 3, Consideremos el PV IF

u(0,x) = f ( x ) , o < x < 1

4.

t iene como representación estocásti ca a

Al respecto, se trata de estimar esperanzas de c i e r t a s V A ' S , por l o que conviene recurr i r a l a Ley de lo s Números Grandes ( L N G ) ; en su f o r - ma fuer te IRuiz Moncayo (1980) I , dicho teorema expresa la consistencia de l a media muestra1 como estimador de 1 a esperanza de una UA, con probabil - i dad 1. En otras palabras, establece la conveniencia de adoptar el siguien - te algoritmo para estimar las esperanzas requeridas para el cálculo aproxi - mado de la solución de un PYF o un PVIF dado:

ALGORITMO

i )F i j a r un punto Q E D y t - > O , si se aplica;

i i )Observar la t rayector ia de 1 a d i fusión que parte de Q

a ) hasta que toque aD por vez erilnera (FVF)

b ) durante t segundos (PV IF)

anotando la posición final P .

¡¡¡)Calcular f ( P ) , donde f e s

a ) l a condición a la f ron tera ( W F )

b ) 1 a condición in i c i a l (PV I F )

iv)Repetir los pasos anteriores N veces, obteniendo los valores es f ( P i ) , i=ly0. . ,N.

v ) Calcular el estimador

5.

A

uN:[f(P1) + ... + f(P,)]/N

Por l a LNG,

h

a ) P ( u N --> u ( Q ) cuando n --> 02 ) = 1 (PYF) ,

A

b ) P ( u N --> u ( t , Q ) cuando n --> m) = 1 (PVIF).

Aunque es razonable en v i s t a de l o ya expuesto, este algoritmo no precisa lo que se entiende por "observar una trayectoria de la difusión", tampoco dice cómo hacerse de una. Lo más recomendable parece ser real i- zar experimentos f i c t i c i o s , con el auxi l io de l a computadora, técnica co - nacida como simulación. En otras palabras, las trayectorias requeridas se generarán medi ante al gori tmos computaci onales adecuados , al desarro- l l o de los cuales dedicamos los siguientes tres capítulos.

En l a monografía [Hernández/Al varez(1985)l hemos presentado l a s técnicas de simulación requeridas para llevar a la práctica dichos algoritmos .

6 .

2. El movimiento browniano

Recordemos 1 a fórmula (1.3) , que resuelve el PVIF Je la sección dnterior, En e l l a f i gu ra una familia de V A ' s , [ X t , t > O sobre un - cierto espacio probabilistic0 que ahí no se precisa , pero que por fuerza está presente. Dicha colección es una instancia par t icular del importan - t e concepto de proceso estocásti co, que permi te ac la rar l a re lac ión en t re l a s Ecuaciones Diferenciales y la Probabilidad.

Brevemente, u n proceso estocásti co ( P E ) sobre un espacio probabi - l í s t i c o dado no es sino una familia cualquiera de VA's sobre dicho espa - cia. A q u í nos interesarán PE's con tantas V A ' s como números reales no negativos, tal como el protagonista de (103): al parámetro t > O se le ident i f icará con e l tiempo. Surgen t a l e s PE's al modelar e l movi - miento de partículas bajo condiciones no determinísticas, como en el s i - guiente ejemplo.

EJEMPLO 1, El efecto Tyndall

Bajo condiciones usuales , el movimiento de una partícula puede ana - l izarse bajo una óptica determinística, es decir , mediante las leyes de l a Mecánica Clásica, Sin embargo, al acercarnos a l o microscópico empiezan a surg i r fenómenos que aparecen como no determinísticos: tal es el caso, por ejemplo , del movimiento de una par t ícula de polvo, ta l como s e ve al iluminarla mediante u n rayo de luz en e l i n t e r io r de u n cuarto obscuro. Las t rayector ias que se observan tienen apariencia errá- t i c a , si bien son continuas , por lo que no parece adecuado modelar e s t e movi - miento mediante el formalismo newtonia- no, En su lugar, construimos u n espa - cio probabilíst i co < ~2. ,a, P > en el

7 .

que R consta de todas las funciones continuas de la va r iab le rea l - t L O, con valores en R 3 , Los eventos son conjuntos de funciones continuas: por ejemplo, el de que " l a s t rayector ia s comiencen de un punto - dado x E R 3 I' es:

{u E R w ( 0 ) = x);

Una convención muy común u t i l i z a l a f u n c i ó n Xt: R -+ IR3 dada por

Xt (u) : = u( t ) , t - > O, con l o que l o s dos eventos anteriores se denotan me -

diante (Xo = x ) y (Xt E A ) , respectivamente.

El efecto Tyndal 1 a que se hace mención en el ejemplo anterior se

conoce también como movimiento browniano, en honor al botánico inglés - Robert Brown, quien lo estudió en un contexto a lgo diferente, en 1828, A l estu - d i a r l a f e r t i l i z a c i ó n de algunas plantas o r i g i n a r i a s de Austra l ia , Brown observó - al microscopio el movimiento de granos de polen suspendidos en agua, fenómeno por de - más e r rá t i co y ajeno a l o conocido hasta - entonces o

Durante buena parte del Siglo X I X fue imposible para los c ient íf i - cos explicar la naturaleza del misterioso movimiento browniano. No fue' s ino hasta 1877 que Delsaux, un f í s ico francés, ofreció una expl icación no v i t a l i s t a , en términos de l a s f ue r za s que sobre l a p a r t í c u l a de polen ejercen 1 as moléculas de agua que 1 a rodean. En 7 L Hernández (1981) 1 nos hemos ocupado del model amiento de este fenómeno desde diferentes ángulos, partiendo de diversas suposiciones acerca de l a con s t i t uc i ón del tiempo y

8.

del espacio. A h í también damos una construcción detallada del modelo de Wiener a l que nos referimos más adel ante.

Einstein, en 1905, dio una teoría cuantitativa del fenómeno (véase [Einstein (1959)] ) , en 1 a que obtuvo l a manera de asignar yrobabi 1 i dades de transición, concepto mediante e l cual se busca predecir el futuro en función de l a información presente; calculamos según

En el tratamiento de Einstein figura la separación del movimiento browniano en t r e s componentes independientes, cada una de las cuales cons

t i tuye un movimiento browniano unidimensio - nal . Físicamente, éste no es o t ra cosa - que l a proyección del movimiento f í s i c o en una de 1 as direcciones coordenadas. Si l o

x, -

t graficamos contra el tiempo, se obtiene - una figura parecida a l a que s e da aquí a 1 a i zqui erda.

Einstein dio l a fórmula para l a densidad de transición en una dimen - sión, a saber

En 1923, Norbert Wiener, matemático de nacionalidad norteamericana, construyó u n espacio probabilístico y un PE definido sobre é1 para los que se cumplen las dos relaciones (1). Dicha construcción constituye u n modelo riguroso sobre el cual se puede basar una t eo r í a del movimiento - browniano que r e su l t e s a t i s f ac to r i a desde el punto de vista matemático; - en honor a su creador, dicho modelo se denomina proceso de Wiener. Breve - mente, el proceso de Wiener es una colección

9.

de VUA's sobre un espacio probabilistic0 4 , ,P>, que sa t i s face

i ) P(Wo = O ) = 1

i i ) Para S < t cualesquiera y

f

con p como en ( I ) ' ' ¶ y

i i i ) Para cualquier elección de n > 1 y de instantes tl,tp, ..., tn

(1 istados en orden ascendente), los incrementos

son V A ' S i ndependi entes.

A p a r t i r de estas tres condiciones, en 1 a teor ía de Wiener del movimiento browniano se prueban dos resultados muy profundos:

a ) P ( t l ->I$ es continua) = 1,

b ) P ( t l ->!E!.+ es diferenciable en algún punto) = O . L.

Así pues, el modelo de Wiener nos da una manera de probabilizar el espacio de las funciones continuas sobre [O,m), según l a cual l a s f u n - cienes que son diferenciables en algún punto caben en un conjunto de probabil idad cero.

Desde el punto de vis ta f ís ico, la cont inuidad de las t rayector ias pone en evidencia el carácter clásico de l a pa r t í cu la browniana en e s t a formulación. Por o t ro l ado , l a no diferenciabil idad de las t rayector ias impide asignar velocidad instantánea a la par t ícu la browniana, l o cual no resu l ta sa t i s fac tor io desde e l punto de v is ta f í s ico . Por esta razón,

10.

los tratamientos físicos acerca del movimiento browniano se basan en - otros modelos, no en el de Wiener-Ei nstei n : por ejemplo , véanse los ar- tículos sobre la ecuación de Langevin en [Wax(1954)] , a s í como [Braun (1984)].

S in embargo, se demuestra en l a l i t e r a t u r a matemática que esos modelos dinámicos del movimiento browniano pueden obtenerse a partir del de Wiener [Nelson(1972)], mediante transformaciones adecuadas de sus t r a - yectori as. Por esa razón el modelo de Wiener continúa en el corazón de la t eor ía matemática del movimiento browniano. Dentro de e s t e e sp í r i t u , en la siguiente sección daremos l a t e o r í a de los procesos de difusión en una dimensión, poniendo énfasis en su construcción a partir del proceso de Wiener.

J

11.

3. Difusiones en una dimensión

En todo lo que sigue <R,a,P> denotará un espacio probabil ístico sobre el cual está definido un proceso de Wiener {Wt,t > O). Para cada t - > O , sea 'la sub- o -álgebra de a generada por los valores pasados y presente de dicho proceso: de manera equivalente, por los incrementos pasados del mismo. En otras palabras ,

t

Observemos que {dt} constituye una fami 1 ia creciente de sub- u - álgebras de& , es decir , una f i l t rac ión , a 1 a cual denotaremos por Q . o

De l a misma manera, usaremos la notación abreviada y, para denotar al PE {xt , t > O). Los PE's en esta sección tienen trayectorias continuas.

Es una consecuencia di recta de las propiedades i i ) y i i i ) del pro - ceso de Wiener (ver la sección 2 ) que

lo cual se abrevia diciendo que el proceso de Wiener es una martingala (MG) [Doob(1953)]. Más aún , esta MG da lugar a una infinidad de otras también asociadas a e s t e proceso, propiedad reproductiva que formará el hilo conductor de cuanto sigue, que mucho debe a l tratamiento dado en [Stroock-Varadhan ( 1 9 7 9 ) j .

Por ejemplo, se prueba también a par t i r de las propiedades i i ) y i i i ) del proceso de Wiener que Wt-t es una MG. De a h í es casi inmedia t o que, s i definimos

2

-

12.

{Mt(f),t 2 O ) es una MG

s i f es cuadrática.

En efecto, basta u n cálculo directo y la consideración de que toda combinación 1 i neal de . MG's es una MG para probar plenamente l a a f i ma - ción anterior.

Ahora b i e n , las funciones dos veces continuamente diferenciables (es decir , "de clase C 2 " ) son aquellas que pueden s e r aproximadas "bien" medi ante u n polinomio de segundo grado, en vi r t u d del Teorema de Taylor. Resulta entonces razonable esperar que ( 2 ) se cumpla para cualquier f E C2. En efecto, se t iene de hecho la siguiente caracterización del proceso de Wiener, para cuya demostración nos referimos a [Stroock-Varadhan op. c i t . ] :

l l U n P E x . es el de Wiener s i y solo si cumple (2), cualquiera - que sea f E c2".

Esta propiedad nos permi t e darnos cuenta del importante papel que en la teoría del proceso de Wiener juega el operador diferencial:

razón por la cual en algunos contextos se llama a este operador el generador del proceso de Wiener. -

Surge de aquí la idea de considerar otros operadores (diferenciales o no) j u n t o con otros PE's, para l o cual introduciremos el siguiente concepto, donde AC denota el espacio de todas las funciones reales, medi- bles y acotadas sobre Rn , normado por el supremo:

13.

DEFINICION 1.

Mt(f):=f(Xt) - f(Xo) - l)f (xS)ds ( 4

Se dice que x. está asociado a A s i M . ( f ) es una MG para cualquier f E C T

En particular, W. está asociado al operador (3) . Asimismo, s i de - finirnos xt:=g(Wt), para g E C 2 , entonces x. está asociado al operador diferencial

En general , se prueba fácilmente que, s i x. está asociado al opera dor di ferenci al

entonces g(x. ) , con g € C 2 , está asociado a

Nos damos cuenta de que 1 a clase de los operadores diferenciales A del tipo ( 5 ) es interesante, pues no se trasciende al transformar sua- vemente procesos asociados a operadores en e l l a . Más aún, s i a ( x ) > O

en (4) y g ‘ no se anula, entonces el operador asociado al proceso trans -

14.

formado por g también tiene coeficiente principal positivo.

DEFINICION 2.

Una difusión es u n PE x. asociado a un operador como (5), con a y b continuos, a positivo.

En otras palabras, el proceso de Wiener es una difusión, igual que sus transformaciones. Queda todavía e l problema de caracterizar de alguna manera operativa 1 a clase de todas 1 as di fusiones y , más importante aún , dar u n procedimiento para construir una difusión dado e l operador corres- pondi ente.

Con ese f in , comencemos por observar que M. ( f ) t iene trayectorias continuas y está adaptado a l a f i l t r a c i ó n a . , que es la generada por los incrementos browni anos. Además, se t rata de una MG de cuadrado integra - ble. En vista del teorema de representación de Clark [Clark(1970)] para cada f E C 2 exis te un proceso medible z. t a l que

i i ) E IzSI2ds < m, i

En i i i ) , l a i n t e g r a l que aparece es una integral estocástica, toma - da en el sentido de Ito. Para es te y .otros conceptos de la teor ía de in- tegración estocástica referimos a los lectores al capítulo 4 de [Wong-Hajek (1985)]

En otras palabras, se cumple que

f t /t f (x,) = f (x,) + ),Af(xS)ds + J,ZdW,

15.

para cualquier f E C2 , relación en l a que aún f a l t a i den t i f i ca r e l i n t e - grando z .

Un resultado muy famoso de Ito afirma que, si x. es una difusión, entonces

por lo cual

f t ít f (x,) = f (x,) +) -Af (x,) ds +]; f ' (X,) Ja(xS)dWS,

O

cualquiera que sea f € C2. En par t icu lar , podemos e l eg i r f ( x ) = x, con lo que ( 6 ) da lugar a

Se acostumbra e sc r ib i r e s t a ir1 tima relación en forma diferencial , como u n problema de valores iniciales .

y entonces se le conoce como una ecuación diferencial estocástica ( E D E ) . A q u i hemos puesto o : = &, notación que conservaremos durante el resto de este trabajo.

De acuerdo con 1 a Regla Diferencial de I to [Wong-Hajek , op. ci t .I , si x. sa t i s face l a E D E ( 7 ) y si f E C2 , entonces f ( x . ) necesariamente sa t i s face (6 ) . Por lo tanto, las difusiones son las soluciones de EDE's - del t ipo de *( 7 ) . En 1 a siguiente sección nos ocupamos brevemente de 1 a construcción de difusiones, a reserva de ocuparnos de dichas cuestiones con mayor de t a l l e en el siguiente capítulo.

-

16.

4, So luc iones probabi l í s t icas de EDP's

Sea x , una difus ión en e l sent ido que definimos en l a s ecc i ón

anter ior, asociada a un operador diferencial A del t ipo (3.5) , con a > O. Sabemos que podemos cons t ru i r 1 as trayector ias de x. transforman - do adecuadamente l a s del proceso de Wiener, según prescribe (3.7). A - continuación veremos cómo expresar las so luc iones de PVF's y PVZF's aso - ciados al operador A como funcionales de la s t rayector ia s de x,

Para comenzar, 1 a ecuaci ó n di ferenci a l estocást i ca

define un algoritmo para construir 1 as trayectorias de x. a par t i r de l a s de W. , a l menos en pr inc ip io . En efecto, podemos cons t ru i r una suce-

s i ó n de aproximaciones {x.(n))> mediante el algoritmo

x o . (n+l): = X +J*b(xs(n))ds +i.a(xs(n))dWs, O O

con x = x, para x E R dado. Se obt iene a s í una sucesión de PE ' s - que, s i b y a son suficientemente suaves converge a un Único 1 ímite; - además, s i l a convergencia es suficientemente fuerte, es automático que e l proceso l ímite sat i s face ( I ) ,

O

Por lo pronto, observamos que, de (2),

donde en los puntos suspensivos no f igura para S > t. P o r l o tan-

to , s e cumple que

17.

Es de esperarse que l a propiedad anterior se preserve bajo la convergencia, si es suficientemente fuerte, De ser as í , se tendrá que

para t > O , h > O a rb i t r a r io s , A un bore1 i ano arb i t ra r io de la rec ta , Si se cumple esta propi edad, se di rá que x. es de Markoff.

-

En el capítulo 4 de [Wong-Haiek op. cit.] se prueba que el método de aproximaciones sucesivas arroja una sucesión convergente de PE's y que todo lo an te r ior se cumple bajo las siguientes condiciones sobre los coefi ci entes de (1 ) :

A ) Condición de crecimiento lineal

B,) Condición de Lipschi t z ,

donde K es una constante positiva adecuada.

Para resumir, bajo las condiciones A ) y B) exis te un fnico PE x. que sa t i s face (1 ), mismo que se obtiene como lími te de l a s aproximaciones sucesivas , en el sentido de que

uniformemente en t si restringimos esta variable a intervalos compactos , digamos [O,T]. A su vez, l a unicidad de la solución así obtenida debe entenderse en el sentido de que, si y. e s o t r a que resul ta de aplicar el algoritmo ( Z ) , entonces

18.

para cualquier T > O .

Más aún, l a convergencia ( 4 ) nos arro ja u n límite x. que es de Markoff, además de tenr trayectorias continuas.

La propiedad de Markoff es muy importante y , en términos informa- l e s , s i g n i f i c a l o siguiente: por lo que se re f iere a predecir el comporta - miento futuro de la t rayector ia , no hay más información en e l pasado que l a contenida en el presente de l a misma. En otras palabras, el presente xt puede servir como condición inic ia l para predecir el futuro,

{xs, S ' t}"

Por lo que se ha dicho, 1 a propiedad de Markoff genera un comporta - miento en el proceso que recuerda el que se observa en las ecuaci ones di - ferenciales. A continuación haremos aún más explícita esta conexión, para l o cual supondremos que se cumplen l as condiciones A) y B ) , de mane - ra que exis te una úni ca difusión asociada al operador diferencial

1 d 2 d dx A:= a(x) - + b(x) - ,

dx

l a cual tiene trayectorias continuas y es de Markoff además de cumplir l a condición in ic ia l x = x. Para hacer más explícita esta última propiedad , denotaremos 1 a correspondiente probabi 1 i dad mediante el símbolo - Px y a l a operación de tomar esperanza por Ex.

O

Bajo estas condiciones, efectivamente tenemos una MG M. ( f ) sobre <R,&,P >, de t a l manera que

X

r.

f (x.) = f ( x ) + J, A f (xs) ds + M. (f)

19.

además de que Mo(f) = O, cualquiera que sea f E C 2 .

En lo que sigue utilizaremos el concepto de tiempo de paro (TP) - asociado a la f i l t rac ión a, Se entiende como t a l a toda V A T que toma valores en O , y , además.

Por el Teorema de Muestre0 Opcional para MG's [ Doob op. c i t j , cualquiera que sea el TP T debe c'umpl i r s e que

En consecuenci a

Exf(xT) = f ( x ) + Af(xs)ds, l: re1 aci ón que s e conoce como fórmul a de Dynkin.

Sea ahora u solución del ?VF

Au = O en (a.b); u(a) = A,u(b) = B ( 7 1

Definamos el tiempo aleatorio T:Q -f [O,-lque marca el instante en que x. sale de dicho intervalo por primera vez, es decir

T(w) := in f { t > O:xt(w) E! (a ,b)) .

Es u n resul tad0 conocido que T es u n TP en el sentido de (5).

Además , es claro que

Au(xs(w)) = O , O < S < T(w) .

A s i xT = a

B si xT = b u(xT) = {

20.

por lo que

Exu(xT) = APx(xT = a ) + BPx(xT = b )

De l a fórmula de Dynkin resulta entonces en forma inmediata que

u(x) = APx(xT = a ) + BPx(xT = b ) , ( 8 ) I

resultado que explica la solución'que dimos del ejemplo 1.1.

La correspondiente formulación para d i fusiones mu1 tidimensionales nos proporciona l a fórmula estocástica para la solución de PvF's como

1 n n

i , j=l j i=l 2 a i j ( x ) axiax a 2 u + C bi (x) = O en D , ax

donde D es u n dominio en Rn y f e s una función medible y acotada en D. Directamente de l a fórmula de D Y n k i n ,

donde T denota ahora el instante en que x. abandona D por primera vez.

Definamos ahora una familia

Vt,t > 01

de operadores lineales sobre AC mediante

21.

P t f (x) := E x f (Xt) (9 1

para f E AC , x E R. Si f está en el dominio de A , de (1) tenemos que

t

O Ptf(x) = f ( x ) +j PtAf (xS)ds.

En e fecto , M. ( f ) e s una MG y entonces su esperanza es constante, es decir, vale cero en es te caso ya que Mo(f) = O . En particular,

P o f = f ,

además de que

- P f ( x ) = PtAf (X) d d t t

para cada f E DA, x E IR . De aquí es claro que

para f y x adecuados , relación que nos da una interpretación del opera - dor A .

Más aún, l a propiedad de Markoff es equivalente a l a condición

Px(Xt+s E dy)P (X E A ) . Y S

relación que se conoce como de Chapman-Kolmogoroff [véase Doob (1953)l . A su vez, ChK se puede escribi r en l a forma de

Pt+s = P t P S (13)"

llamada l a propiedad de semigrupo de { P t l .

22.

Por 1 a propiedad de semigrupo, es claro que los operadores de l a fami 1 i a { Pt} conmutan. De ahí es inmediato que

P h f -f Ph ( Pt f )- Ptf P t ( T ) = h

y , tomando el 1 ími te cuando h -+ O , nos convencemos de 1 o sigui ente: P t f E DA para cualquier f E AC y , además,

APt f = PtAf. \

En par t icular , de (11) obtenemos que

- P f ( X ) = APtf(x) d d t t

Las relaciones (11) y (14) se conocen como ecuaciones de Kolmogoroff, hacia adel ante y hacia atrás, respectivamente y forman l a base del enfoque tradicional al estudio de los procesos de difusión. Para nosotros, la ecua - ción (14), j u n t o con (9), nos indica que l a función

u ( t , x ) : = P t f (x )

resuelve un PVI. Más especificamente, la solución del PVI

" - A U en (O,..) X D a t u(0.x) = f ( x ) , x E D

t iene como representación estocástica a

u ( t , x ) = E x f ( X t ) .

Con e s to queda establecido sobre bases firmes el ejemplo 1 de l a

23.

sección 3 y terminamos nuestra breve presentación del tema. Recomenda- mos a los lectores interesados consultar el capítulo 6 de [Friedman (1975)l para ampliar sus conocimientos sobre el tema,

1.

C A P I T U L O

Introducción

Estamos interesados en la so luc ión de l a EDE(1.3.7), que por como-

didad volvemos a escribir:

dxt = b(xt)dt .f o(xt)dWt; x. dado

Desgraciadamente , só l o en muy contados casos es fact ible resolver-

l a analíticamente, y encon t ra r a s í l a d i f u s i ón x. que l a s a t i s f a ce . En

l a mayoria de los casos deberemos contentarnos con encontrar aproximaciones numéricas de x. Más precisamente, s i contamos con algún procedimien - to numérico para generar trayectorias aproximadas de movimiento browniano W. (ver el capítulo 3) aunado a un método de integración numérica, podremos generar trayectorias aproximadas de 1 a di fusión x. Así , para cada - trayectoria aproximada de W. obtendremos una trayectoria asociada de x.

Repitiendo el procedimiento un c ie r to número de veces, podremos ob - tener una muestra de trayectorias aproximadas de l a d i f u s i ó n x. en un intervalo de tiempo, digamos [O,T] , que puede ser aprovechada de var ia s ma- neras. Por ejemplo, para algún instante de interés LE[O,TJ tendremos - una muestra de la va r iab le a leator ia x t , que nos permi ti rá e s t imar su s propiedades distribucionales, O bien, podemos estudiar tiempos de paro, por ejemplo el tiempo T que tarda el proceso en sal ir del intervalo (a,b) - partiendo de x. E (a,b). O alternativamente, .bajo l a s mismas condiciones podemos es t imar l a p robab i l idad de que alcance el valor a antes que

e l b.

Tendremos dos fuentes de error a l generar las trayector ias aproximadas de l a d i f u s i ó n x . : l a primera resulta de l a d i s c r e t i z a c i ón del - tiempo y el método de integración numérica que usemos; l a segunda de usar una aproximación del movimiento browni ano en 1 ugar del proceso real W . En las secciones 2 y 3 analizaremos por separado estas dos fuentes de - error, y en l a s e c c i ón 4 mostraremos algunos resultados computacionales.

25.

2. Efecto de la discretización del tiempo.

En esta sección supondremos que contamos con trayectorias exactas de movimiento browni ano W . y no sólo con aproximaciones.

Aprovechando el carácter markoffiano de l a d i f u s i ó n x. buscaremos soluciones aproximadas x. de ( l o l ) d e f i n i d a s recursivamente en u n conjun to discreto de tiempos:

- -

o = t o < t < . o . . . 1 = "

empezando en yo = x. y mediante algún método de integración numérica calculamos las siguientes xi como aproximaciones de x t i 9 i = 1 , 2 9 . 0 . 9 n . En los puntos intermedios usaremos interpolación l ineal:

-

obteniendo así trayectorias continuas xt en [ 0,TI. Por simp1 i cidad, en l o que sigue supondremos que tenemos incrementos de tiempo iguales: titi - ti = h; i = 1 , 2 , . 0 . , n .

-

Un esquema senc i l lo de integración numérica es el análogo estocás- t i c0 del método de Euler:

donde AWi representa el incremento de W t en el intervalo [ti ,titi].

Maruyama (1955) demostró que e s t e método converge uniformemente en media cuadrática a x t 9 o sea que:

26.

Podemos obtener mayor precisión con e l método de integración de Heun, que resul ta de uti l izar la regla trapezoidal en vez de l a escalona - da para aproximar las integrales sobre [ti, ti+l]:

con el predictor de Euler:

McShane (1974) demostró que e s t e método también converge uniformemente en media cuadrática, en el mismo sentido ( 2 ) , pero ya no a l a solu- ción de (1.1). En efecto, McShane prueba que e l proceso 1 ími t e x. s a t i s face :

-

Nótese que (1.1) y ( 4 ) sólo difieren por el factor que multiplica a d t , por lo que a l u t i 1 i zar e l método de Heun (3) debemos tener cuidado de introducir u n término de corrección en b (x t ) . De es ta manera, el mé- todo de Heun (3) da u n resul tad0 equivalente al método de Euler (1).

Observando que los métodos de E u l e r y He& son ambos del t ipo Runge-Kutta, Rümelin (1952) estudió los esquemas de este t ipo para el caso estocásti co, abarcando a s í una c lase muy amplia de métodos numéri- cos. En l a notación de RUmelin, los métodos de Runge-Kutta de orden m + l son como sigue:

m m - 'i+l -

- - + C p . K . h + C q . G . AWi j =O J J j =O J J (5) '

27.

donde :

KO = b ( x i ) Go = a(';;i)

x i ( ' )= Xi + Bl0 KO h + yl0 Go AWi

K1 = b(xi (1)

. . . . . . . . . . .

con la restr icc ión:

m m c p j = c q j = l

j =O j = O (5) I l '

Basándose en e l resultado de Maruyama (1955), Rtimelin (1982) demos - trÓ que cada uno de estos métodos Runge-Kutta de orden mtl también converge uniformemente en media cuadrática, en el mismo sentido (Z ) , pero a l a solución de:

dx t = [b(xt) + A&(, dx t ) a(xt)]dt + a(xt)dWt (6) '

con el factor de corrección X dado por:

x =(" si m=O

(6 ) "

Lo Único que se requiere para que la solución aproximada x. que - arroja el esquema (5) converja a la solución de l a ecuación ( 6 ) es que:

sean acotadas. Nótese que este requerimiento no representa rljnguna res- tricción adicional cuando trabajamos en una computadora, ya que cualquier número generado en el 1 a e s t á acotado por su propi a capacidad [ Rümel i n (1982)] .

Para obtener una solución de l a EDE (1.1) con el método Runge-Kutta ( 5 ) , sólo debemos sustituir el término b ( x t ) por el término corregido:

en (5), ya que entonces (6 ) equivale a (1.1).

RUmelin (1982) también investigó el orden de l a convergencia, y l legó a la conclusión de que e l esquema de integración numérica (5 ) a l - canza el mejor orden posible,

Para precisar, diremos que se t iene u n orden de convergencia O ( h r ) si

E (Ft - x ~ ) ~ = O ( h r ) para t E [O,T]

donde

29.

lim sup ( l / h r ) O ( h r ) < c o o

h + O

El orden de convergencia que puede ser alcanzado mediante integra - ción numérica depende de l a función:

Si S # O , el mejor orden,posible es O(h2) y l o alcanza el méto- do de Heun (3), uno de la famil ia (5), para e l cual X = 1/2.

Si S = O , el mejor orden posible es al menos O ( h 3 ) . Este orden l o alcanza el método c lás ico de Runge-Kutta de cuarto orden (m = 3 ) , con tenido en (5), para el cual X = 1/2.

-

Para l legar a estos resultados, se supone que b y CJ y sus derivadas parciales hasta de cuarto orden son acotadas y continuas.

Estos resultados pueden extenderse al caso en que b y o no sólo dependen de xt , s ino también de t :

b = b(t ,xt ) u = u ( t ) 'Xt

También pueden extenderse al caso en que x. y W. sean mu1 tidimen - sionales. Al lector interesado en estas extensiones se le recomienda con- su l tar [Rümelin (1982)].

30.

3. Efecto de l a aproximación del movimiento Browniano.

En esta sección anal izaremos el error que se produce en 1 as trayec - torias aproximadas de l a di fusión x. cuando usamos 1 as trayectorias aproxi - madas del movimiento browniano W. , siguiendo el t raba jo de Janssen (1984), al que remitimos a los lectores interesados en los detal les.

En l a EDE ( 1 . 1 ) supongamos que los coeficientes b y CJ sa t i s fa - cen una condición de Lipschitz uniforme y que CJ es acotada,

Por los resultados de l a sección anterior sabemos que 1 a solución de l a ecuación en diferencias:

- x. = x ; n = O , l , ..., N-1

Converge uniformemente en medi a cuadráti ca a xt. Al implantar en l a computadora el esquema en diferencias ( 1 ) , obtenemos trayectorias aproximadas del proceso de Wiener W. Lo que se ha hecho entonces es - reemplazar W , o t r a variable aleatoria W resultando en consecuen- c ia u n nuevo esquema en diferencias:

n ,m

X 0 ,m = x ; n = O , l , . .. , N-1 (2 1

s i ~ w o Y m ~ . * u w 1 4 ~ ~ , m ) 3 (Awl ,. . . , AM, d ) donde D denota convergen-

cia débil , entonces X -P Xn. En efecto, es fáci l ver que D n ,m,,,

X - -

n ,m - Wo,m.Wlym,. . . "N-1 ,m) donde $, es una función -

31.

continua; de aquí la afirmación. Entonces estamos interesados en obtener estimaciones para la rapidez de dicha convergencia, La respuesta la dare - mes en términos de la distancia Lipschitz cuya definición se da a continuaci ó n o

tri bución de Xn+l (resp. x ). Por 1-1 (resp. pm) n + l .m w 1. n,m

Siguiendo a Janssen L13841, denotaremos por LI (resp. l a dis- denotaremos l a dis-

Definamos , para c&

L i pschi t z como:

i nimos l a correspond iente dis tancia de

son distribuciones

de probabilidad en Rj . Observemos que esta distancia está acotada ya que g , IT^, n 2 , lo están. Además, es natural introducir esta distancia ya que

estamos interesados en 1 a rapidez de 1 a convergencia de X D '3' m+. m

+ xn.

Para simplificar la notación pondremos:

Por ú1 timo, L denotará la constante de Lipschitz para b y u y M una cota superior para o.

Tenemos así el siguiente resul tad0 cuya demostración se encuentra en Janssen [1984]:

32.

Teorema 1. Sean W var iables a leator ias independientes tales j , m

que E 1 Wj - < h 5 1 o Entonces, existe una constante C1 que depende

só l o de L, ta l que para toda r > O:

T d(p,pm) 5 4r-2 exp(clhi) + r(r+l) 5 (y,Fm).

Denotemos por Vh l a d i s t r ibuc ión Gauss iana de media cero y var ian - za h., puesto que la s va r iab les a l eator i a s AWi son independientes y - gaussianas podemos expresar a ji como un producto de medidas unidimensio- nales:

- p = vh + O 0 0 + Vh. Supongamos que s e descompone de 1 a misma m

manera :

- - %,m - - + .. . + En,,,. Obtenemos entonces que:

n d(L, Fm) 5 c d(Vh, pi ,m). Si suponemos además que l a s W. son iden

ticamente distribuidas, entonces 1 as pi ,m son igua les y por l o tanto

- i =O 1 ,m

-

En el contexto del teorema 1, es to s i gn i f i ca que para garanti zar un e r ro r de di scretización del tiempo de tamaño O (h) , necesitamos una - aproximación de V,, de tamaño O(h2), es dec i r , d(Vh, = O(h2).

-

Solo nos hace fa l ta e s t imar l a d i s tanc ia ent re Xt y

XN,m(t):=Xn ,m (t) s i t E [tn,tntl).

E l s igu iente teorema nos proporciona tal estimación.

33.

Teorema 2. Sea t E tn+2 ) . Denotemos por px (res p. t

pX ( t ) ) las dis tr ibuciones de xt(resp. x ( t ) ) . Entonces:

N ,m N ,m

Recordemos que el teorema 1 nos da una estimación de d(p,pm) y

a s í obtenemos una estimación del error que cometemos al aproximar la di - fusión xt por medio del esquema en diferencias (2).

34.

4. Resultados computacionales.

Se implementó un algori tmo computacional para encontrar trayecto- r i a s aproximadas de la difusión x. que s a t i s f a c e l a EDE (1.1). La gene - ración de trayectorias aproximadas del movimiento browniano W. se hizo - con las técnicas que se describen en el capí tulo 3, basadas en caminatas aleatori as ( C A ) y procesos de transporte uniforme (PTU). El método de - integración numérica de Runge-Kutta usado es el desarrollado por E. Fehlberg (1970) e implementado con la ru t ina RKF45 por L,F. Shampine y H. A. Watts (1974) [Forsythe e t ' a l . (1977)].

E l algoritmo se probó para el caso particular en que:

a ( x ) = J x ( 1 - X ) / n

donde S y N son constantes. En e s t e caso se tomaron:

S = 0.1, N = 10

Así pues a 1 a EDE (1.1) se convierte en:

dx = sx ( l - x t ) d t + 4~ (l-xt)/2N t t t dWt; x0 dado

Se recomienda ver el capítulo 5 para una interpretación biológica de ( 2 ) .

En 1 as gráficas i a 3 se muestran algunas de 1 as trayectorias aproximadas de x. que se obtienen con este algoritmo, que Se pwx?nta en e l Apéndice B.

Para cada uno de los cinco puntos iniciales x. = 0.1,0.3, ... ,0.9

se generaron cien trayectori as aproximadas de x. usando PTU, y cien más usando C A , en el intervalo O S t 5 50. En l as g ráf icas 4 a 8 se - muestran los h i stogramas con 1 as frecuencias observadas de las var iables

35.

a lea tor ias x t , para t = 5, 10, 1 5 , 20, 35, 50. Nótese que en todos

los casos obtenemos resultados similares usando PTU o CA,

Estas frecuencias observadas tienen una interpretación bien preci - sa. Son estimadores de las probabilidades de transición P(xt E (a ,b) I xo) , que en conjunto caracterizan completamente al proceso x.

Con esta herramienta cornputacional también s e pueden estudiar tiem - pos de paro , como por ejemplo:

con :

a = 0.02 y b = 0.98 (3)"

Para cada uno de los nueve puntos i n i c i a l e s x. = 0.1,0.2, ..., 0.9

se generaron cien trayectorias usando PTU, y cien más usando CA, parán- dose el proceso en el tiempo de paro definido en (3) . En 1 a gráfica 9 se muestra l a esperanza estimada de T , como función del punto in ic ia l xo.

En 1 as gráf icas 10 a 12 se muestran los h i stogramas con l a s frecuen - cias observadas de T o Nótese l a s imi l i t ud de los resultados obtenidos usando PTU o CA.

En el capítulo 1 se mostró que l a EDE (1.3.7) está relacionada con la solución del PVF (1.4.7). Así como en nuestro caso particular la EDE (1.3,7) se convirt ió en ( 2 ) por la sust i tución ( l ) , e l PVF (1.4.7) se convierte en:

1 x (1 -x ) d2u - ( x ) + sx(1-x) ay (x) = o du 2 dX2

u(a) = A , u(b) = B

36.

que se puede simplificar:

Fácilmente se puede ver i f icar que su solución está dada por:

donde :

g(x) := exp(-4Nsx) (5)”

Ahora, con T definido como en (3), l a ecuación (1.4.8) nos conec- t a l a EDE (2) con el PVF ( 4 ) :

Comparando ( 5 ) y (6 ) tenemos :

Esta re1 ación es de gran importancia para nosotros, porque nos per mi te verif icar el algoritmo usado para la solución de (2) o Para tal efec- t o podemos recurrir a los estimadores de Px(xT = a ) que nos proporcionan las trayectorias generadas con anterioridad. En 1 a gráfica 13 se mues - tran estos estimadores (con sus respectivos intervalos de confianza al 95%) superpuestos al valor exacto calculado por medio de l a ecuación ( 7 )

-

37.

Nótese que casi todos los intervalos de confianza al 95% para l a estimación de Px(xT = a ) engloban el valor calculado por medio de l a

ecuación ( 7 ) , lo cual nos indica que hemos obtenido u n resul tado sat is- factorio.

G R A F I C A 2

L

i i

41.

U 5 1

0.0 1 6 2

U

PTU t=5

GRAFI CA 4-a

Frecuencias de Xt, para X. = O. 1

0.2 0.4 0 - 6 0.8

PTU t=10

0.0 0.2 0.4 0 . 6 o . 8

6 5

U

o . o

PTU t= 15

'O . 2 0 . 4 0 . 6 O . 8 1 . o

ij

O . I 0 . 0 I 1

o 1 . o

CA t=5

12 31% 0.2 0 . 4 0.6

ü

1

CA t=10

1 o . 2 0 - 4 2 0 . 6 0 . 8 1. o

CA t= 15

"I 1 0

u . 2 0 . 4 0 . 6 0 . 8 1 . 0

42.

6 6

-1 GRAFICA 4-b

Frecuencias de X t , para X. = 0.1

P TU t=20

13.

n 3 3 3 3 4 3

2 l l l l l ” - - o. o 0 . 2 0.4 0 . 6 0 . 8 1 . 0

6 7

PTU t=35

4 4

0 . 0 0.2 0 . 4 0 . 6 0 . 8

u PTU t=50

‘1 1.0

0 . 0 0.2 0.4 0 . 6 0 . 8 1.0

o . o

7 1

CA t=20

O

CA t=35

1 6

o. 2 0.4 0 . 6 0 . 8 1 . o 1 I 1

j 0 . 8 1 . o

o . o

CA t=50

1 1 1 1 1 1 1

o. 2 0.4 0 . 6 0 . 8 I 1.0

43.

PTU t = 5

GRAFIcA 5-a

Frecuencia de X t , para X. = 0.3

0 . 0 0.2 0.4 0 . 6 0 . 8 1 . 0

PTU t=10

CA t = 5

2 0

0 . 0 0 . 2 0 .4 0 . 6 0.8

2 1

n CA t=10

1 5 - 16

12 11 -

10 10- 1 2 1 1 “ 1 0 9 8 9 9 - 4

6 6 6 2

c

- o. o 0.2 0.4 0.6 o . e 1.0 o .o 0.2 0.4 0.6 0.8 1.0

n 3 0

PTU 20 - t=15

6 - 8 -

0 . 0 0 . 2 O . * 0 . 6 0 . 8 l.( 3

n. CA t= 15 n n

11 1 2 1 ._

7 9

- 5 6 -

- Lt - 0.0 0.2 0.4 0.6 0.8 1.0

44.

GRAFICA 5-b

Frecuencias de Xt, para X. = O . 3

3 9

ri PTU t=20

o . o 0 . 2 O . + 0 . 6 0 . 8 1 . o

0 . 0 I ü

2 6

o . o

PTU t=35

0 . 2 0 . 4 O .6 O .8 1 . 0 o . o

2 5 PTU t=50

3 1 ' 5 I 5

0 . 2 O . 4 0 . 6 0.8 I

4 9

CA t= 35

0.2 0.4 O . 6 0 . 8 1 . 0

6 3

CA t=50

1

o . o 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1.0

. o

45.

GRAFICA 6-a

Frecuencias de Xt, para X. = 0.5

PTU t=5

1 9

0 . 0 0.2 0 . 4 0.6 0 . 8 1. o

P TU t=10

1 7 i 0 . 0 0.2 0 . 4 0 . 6 0 . 8 1 I

CA t=5

n 20

I 1s

I 0 . 0 0.2 0 . 4 O .6 O . 8 1.0

CA t= 10

i K i , , m 4

2 2

0 . 0 0.2 0 . 4 O . 6 O .8

CA t=15

I 1.0

4 8

U - 1 2

7 5 4

6 6

2 2 I I '

0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1.0

46.

GRAFI CA 6-b

Frecuencias de X t , para X. = 0.5

PTU t = 2 0

8

o . o 0 . 2 0.4 0 . 6 0 . 8

PTU t = 3 5

0 . 0 0 . 2 O . ' + 0 . 6 0 . 8

PTU t = 5 0

] I 1 . o

c

- 1 . o

4

CA t = 2 0

1 1 . I I 5

1 1

o . o 0 . 2 0.1, 0 . 6 0 . 8

CA t=35

, l , l ,

0 . 0 0 . 2 0.4 0 . 6 0 . 8

CA t=50

1 3

I

47.

GRAFICA 7-a

Frecuencias de Xt, para X. = O. 7

3 5

PTU t=5

- 1 8 1

u. 2 0 . 4 0 . 6 0 . 8

PTU t= 10

r"

o . 2 0 . 4 0 . 6 0 . 8

1 1 . u

PTU t=15

* 2 2 2 3

I l l r - o. 2 0 . 4 0 . 6 0 . 8

t . 1 . u

CA t=5

0 . 2 0 . 4 U . 6 0 . 8

CA t= 10

0 . 0 0 . 2 O . 4 0 . 6 9 . 8

I i

CA t= 15

6 8

48.

PTU t=20

GRAFICA 7-b

Frecuencias de Xt, pa ra X. = 0.7

7 5

3 3 2 2 2

0 . 0 0.2 0 . 4 0 . 6 0 . 8 1 .o

a 7

PTU t=35

CA t=20

f r 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . O

CA t= 35

0 . 0 0 . 2 0.4 0 . 6 0 . 8 1 - 0 O . O 0 . 2 0 . 4 O . 6 U - 8 1. o

PTU t=50

CA t=50

49.

GRAFICA 8-a

F r e c u e n c i a de Xt, p a r a X. = 0.9

PT U t= 5

P TU t= 10

7 1

Ir

8

3

- 0 . 6 0 . 8 1 . o

8 2

U

9 1

CA t= 5

- 9

U - 6 0 . 8

CA t=10

8 0

P

O - 4 0 . 6 0 . 8 1 . 0

50.

GRAFI CA 8-b

Frecuencia de Xt , para X. = O. 9

P TU t=20

0 . 6 0 . 8

9 7 3 1 . o

9 9

P TU

t=3*

PTU t=50

0 - . 8 1 . o

9 9

1

CA t=20

1 2 1 1 2 t i !

0 . 0 o ,2 o .4 0 . 6 0 . 8

CA t= 35

o . o 0 - 2 0.4 0 . 6 0.8

9 3

c

1 .u

9 7

- 1 . o

CA t= 50

2 i A

o D O 0 . 2 0.4 0 . 6 0 . 8 1 - 0

CA t= 50 N

o D O 0 . 2 0.4 0 . 6 0 . 8 1 - 0

GRAFICA 9

51.

A N - r )

3 0 - ,

2 0 -

1 0 -

o 1 , I I t> a O. 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 0 . 7 0 . 8 0.9

xO

D Ud

. 3 0 -

2 0 -

1 o "

O ' - I I

a 0 . 1 0.2 0 . 3 0.4 0 . 5 0 . 6 0 . 7 0 . 8 o . 9

52.

GRAFICA 10

Frecuencias de T = i n f {t 0:X E f ( a , b ) }

L a X,=b; b=O. 98

O 20 4 0 6 0 7 0

3 9

O 2 0 40 6 0 m

if I t 7 PTU

I I I 1 6

x. = 0.3

CA x* = 0.2

3 1

Ly- . ...... -x= 1 1 o 1 0 3 0 5 0 7 0 8

n

53.

G R A F I CA 11

Frecuencias de T:= i n f {t. > 0:X d ( a , b ) ) - T

m x = 0.02 T

T 1 PTU

u X = 0.98

I L x. = 0.4

n

O 2 0 4 0 6 0 7 O

3 7

o I l l 2 0 4 0 6 0 co

a = 0.02

b = 0.98

,

O 2 0

CA x. = 0.5

3 5

O 2 0 4 0 6 0 7 0

54.

G R A F I C A 12

Frecuencias de T : = inf{t - > O: X T L(a,b))

m x = a; a=0.02

0 X = b; b=0.98

P TU X. = 0.7

O 2 0 6 0 7 0

PTU x. = 0.8

3 2

8 4

P TU x. = 0.9

U

O 2 0 4 0 6 0 7 0

5 2

CA X. = 0.7

7

1%

O 2 0 40 6 0

O i 2

86

CA X. = 0.8

& 3 0 5 0

CA x. = 0.9

7

O 2 0 ‘10

55.

GRAFICA 13

Px(XT= a ) ; a = 0.02, b = 0.98

PTU

a 0.1 0 . 2 0 . 3 0 . 4 0.5 0 . 6 0 6 7 0 . 8 O 1 9 b X

CA

C A P I T U L O I11

3. l . Introducción

En el capítulo anterior se resolvió la E D E .

dxt = b ( x t ) d t + (x,) dW ; x. dado (1)

u t i 1 izando un método numérico, y fue necesario uti 1 izar t rayector ias de movimiento browniano (b1B). En este capítulo veremos l a simulación de es tas trayectori as.

-

Según ya señal amos en la sección 1.2, casi todas 1 as trayectorias W(u) de MB son continuas pero no diferenciables. [Loeve,M.(1977)]. Esto quiere decir que para todo instante t > O , al observar el comportami ento de cada trayectoria W( 0 ) en cual quier vecindad ( t - E , t + E ) de t , dicha - trayectoria es bastante irregular. A s í pues, si pensamos simular trayec- to r ias de MI3 tendremos que recur r i r a alguna familia de procesos es t o - cásticos que se aproxime al MB y sea tal que sus elementos (procesos) tengan trayectorias que puedan generarse fácilmente.

Si vamos a simular trayectorias de MB a través de otro proceso, es necesario saber qué propiedades se deben sa t i s facer ;

Estas propiedades son las que caracterizan al M6 de acuerdo a l a siguiente

DEFINICION 1.

E l proceso estocástico W = { N (u) ; t > O , u E a} es u n MB s i sat isface:

a ) W t iene incrementos independientes; e s to e s , s i n es cualquier entero positivo y s i O = t < tl < . .. < tn < m, entonces los n incre-

rnentos O

5 7 .

{V(t .) - W(tj-l), 1 < j < n } J

son variables aleatorias independientes.

b ) La distribución de W(t) - W(s) e s normal con media cero varian - za ( t -s)a2, donde O < S < t <a y no depende ni de S n i de t.

c ) PtW ( O ) = O ) = 1

En las siguientes dos secciones haremos ver que podemos aproximar- nos al MB a través de caminatas a leator ias y también a través de Proce- sos de Transporte Uniforme.

3.2 Aproximación de movimiento browniano mediante caminatas aleatori as

Construiremos u n proceso senci l lo cuyas trayectorias sean fáciles de calcular y que se aproxime a un MB. Tal proceso será una caminata a leator ia .

Sea el proceso { X t I t - > O } t a l que X. = O y que a incrementos de

tiempo T t iene sal tos independientes, dados por una variable aleatoria Bernoulli Z. Este proceso X es ta rá dado por:

donde Zo = O

y Z1, Z 2 , ... son copias independientes de Z , cuya distribución está da -

da por:

58.

A l proceso anterior {Xt I t > O ) s e l e conoce como caminata aleato-

r i a (CA)

Nos interesa ver bajo que condiciones el proceso X anter ior converge a un MB, para l o cual hay que trabajar con 1 a d i s t r ibuc ión de Xt.

Cada sal to Z es una v.a. con l a d i s t r i b u c i ó n dada en (Z ) , a s í

que su funcidn generatr iz de momentos (fgm) es:

mz(0) = E(eez) = peeh + qe -0 h

En el instante t se habrán dado n = [t/TI s a l t o s , a s í que e l - desplazamiento total será:

n

i =1 Xt = c zi

y como Z1, Z2,. . .Z son va. i .e.d. tendremos que l a fgm. de Xt será: n

De l a fgm. m ( O ) podemos obtener los dos primeros monentos de Xt: Xt

V(Xt) = m" (O) - (m;( ( O ) ) 2 = 4npqh2 x+ + L L

Nos interesa que cuando T -t O y h -+ O, 1 a CA Xt t ienda a un MB. Esta exigencia, en términos de los dos primeros momentos, tanto de l a CA como del MB se traducen en :

E(Xt) + O, V(Xt) -t to2 cuando T + O y h + O

59.

ya que un MB tiene media cero y varianza to2. Si además tomamos n = [ t / T j ,entonces debe tenerse que:

y cuando T -+ O y h + O

Se tiene el primer límite si tomamos p = q = 1/2, ya que a s í E ( X t ) = O y V(Xt) = ( t / T ) h 2

de aquí que s i T y h tienden a cero de tal forma que h 2 / T = o2 entonces V ( X t ) = to2 . Por l o tanto, el proceso Xt dado por:

ct 1.1 Xt = c zi

i =O

coincide en los dos primeros momentos con un MB.

Otra coincidencia del proceso definido por (3) con un MB Xt

radica en que es un proceso de incrementos independientes.

En efecto, para instantes cualesquiera tl< . .. < tn los incrementos X( t l ) , X( tpl - X ( t l ) ’ . . .w n ) - X(tn-l)

son sumas parciales de diferentes Zi ‘ S , y como estas son Z independien -

tes , entonces también lo son los incrementos del proceso Xt definido

por (3).

Hasta ahora s6lo f a l t a ver que e l l ímite de Xt cuando T -+ O , ten - g a una distribución Normal , esto es :

60.

D -f N(0,a2t) cuando r + O Xt

En efecto , como podemos probar, 1 a fgrn. de Xt converge hacia

la fgm de l a distribución N(0,r2t), ya que s i n = t/.c entonces h = c~fi

V

cuando T + O , donde Y - N(0,a2t).

La convergencia anterior no es mas que el teorema central del 1 í- mi t e , con lo que hemos probado que el proceso Xt converge en distr i bu- ción a u n MB, con varianza 02t.

La confianza de que las trayectorias de CA se aproximan a trayec - torias de MB no descansa únicamente en l a convergencia en distribución que hemos visto, s ino además en el hecho, demostrado por [Knight, F.(1981)1 de que:

Con probabilidad 1 , l a sucesión de CA converge uniformemente a MB .

3.3. Simulación de movimiento browniano por caminatas aleatori as

Si nos interesa muestrear una trayectoria de MB simulado por CA a intervalos de tiempo constantes DT, solamente habrá que generar una trayectoria con punto de partida en cero y ver que valor toma esta en los instantes DT, Z D T , 3DT, ... . De acuer% con la sección 2, los saltos de l a caminata deberán s e r de u n valor A"pequeño", comparado con el intervalo de muestre0 DT, lo cual es equivalente a que l a caminata tenga un número n "grande" de sal tos en el intervalo de tiempo DT, donde n = [DT/A~. Diremos que 1 a CA generada es MB, cuando satisfaga las pro piedades ( a ) , ( b j y (c) de la definicidn 1.1. Ahora bien las propiedades

61.

(a) y (c) se cumplen siempre, ya que l a caminata generada se i.nicia en ce ro y además los saltos de ésta son independientes. Asi que solamente habrá que ver para que valor de n se cumple ( b ) , e s t o e s , para que valor de n , los incrementos

x n , xpn - x n , X3n - xpn 9

%on I' norma 1 es.

Podemos obtener una cota in ic ia l de n utilizando el siguiente teorema de Berry- Esséen [Fe1 1 e r , W. (1968)l:

TEOREMA 1.

Sean Z va. independ que :

ientes, con una distr ibución común F t a l

Sea Fn 1 a distribución de 1 a suma normal i zada

(Z1 + z2 + . . . . zn)/ufi

Entonces , para toda n :

donde rl es la distribución normal estándar.

En nuestro caso, 1 as variables Z son ta les que

p{z , = +h) = P{Z, = -hl = 1f2 , por lo cual

62.

E(Zi) = h2(1/2) + h2(1/2) = h 2 = .a2

E( I Z K I 3 ) = h 3 ( l / 2 ) + h 3 ( l / 2 ) = h3 = p

Por o t r o lado, a los incrementos x j - xj-l = z1 + z2 + . .. + 'n

l e s corresponde l a función de distribución F n , as í que s i estamos dispues -

tos a to lerar una distancia d entre l a función de distribución empírica Fn y la normal estandar,entonces

(3p)/(a3fi) < d.

De aquí tenemos que e l nGmero de sal tos de l a caminata en cada incremento del MB simulado deberá ser :

n > ( 3 p / c ~ ~ d ) ~ = (3/d)2

En la siguiente tabla tenemos los valores de n correspondientes a diferentes valores de d:

d 0.2 O . 15 o. 10 0.05 0.01

n 225 40 O 900 3600 90000

Debido a l a general i d a d del teorema anterior, e l valor obtenido de n resulta innecesariamente grande. Por e l l o , para ver que valor de n es suficientemente grande seguiremos un cr i ter io es tadíst ico que se mencio - na a continuación.

Para u n valor de n dado, se simularon N trayectorias de MB

con M intervalos de muestre0 de n s-a1 tos de l a CA por intervalo, obte -

63.

niendo as? N trayectortas

{M(j,K), j = 1,2 ,... ,M) K = 1 , Z ¶.. .¶N.

Luego se formaron M muestras de tamaño N con los incrementos

{ I(j¶I<): = W(j,K) - W(j-.l,KJ; j=1, ... ,N) y k=l, ..., M

y para cada muestra se probó l a hipstesis nula:

Ho:La muestra { I(j,K); k=l,. .. ,N) proviene de una distribución N(0,DT) contra la hipótesis alternativa:

Hi :Ho es falsa.

La hipótesis de normal idad se probó uti 1 izando 1 as siguientes prue - bas [Seber, G.A.F. (1984)j :

l. Prueba de D' Agosti no 2. Prueba de sesgo de D' Agostino-Pearson 3. Prueba de Kurtosis

Se efectuaron las pruebas para 100 muestras y se registró 1 a propor - ción de rechazos obtenida, 1 a cual se presenta en la tabl a 2 para las 3

' pruebas en el orden correspondiente

-~ ~

.n 10 20 25 30 35 40 45 50 55 60 """""""""""""" """"""""~"""-"--"--"-"------"- """"""""""""""""""""""""""""""""""""" P1 ( n > .43 .20 .13 .o9 .12 .O8 .O3 .O4 .o2 .o9 P2 ( n ) .29 .I9 .16 .O7 .14 .11 .O6 .O6 .o2 .o7

P3 ( n ) .35 .13 15 .O6 .16 .O4 . O4 .O3 04 .o6 """"""""""""."""""""""""""""""""""""""

TABLA #2

En l a tabl a 2 vemos que para u n número n > 40 1 as proporciones

64.

empiezan a osc i lar alrededor de,

Del resultado de las pruebas efectuadas, tenemos que podemos simu - 1 a r MB con CA tomando n > 40 sal tos de l a caminata por cada inter- val o de muestre0 de 1 a trayectoria del MB.

Para l a simulación del MB se u t i l i z ó l a subr1.rtina CAMINA, que ge - nera incrementos independientes de MB para valores deseados de los pará - metros, t a l y como se describe en e l apéndice A .

3.4. Generación de movimiento browniano mediante procesos de trans porte u n i forme

Otra forma de construir un proceso de MB es tomando el 1 ímite de una sucesión de Procesos de Transporte Uniforme (PTU).

Para u n entero positivo n , sea el proceso estocástico v ( t ; n ) una Cadena de Markoff a : Tiempo Continuo (CMTC), con espacio de estados E = i - n , + n ) y generador infinitesimal :

A partir de v ( t , n ) obtenemos e l proceso integrado

W ( t , n ) =I v(t,n)ds 't

O

que se conoce como PTU.

En el capítulo 4 pueden encontrarse mayores detalles acerca de los conceptos recién mencionados.

La jus-t i f icación de esta construcción de MB, está dada en e l s i - guiente teorema, demostrado en:

[Griego et al (1971)l:

TEOREMA 1.

65.

Sea {v(t,n)) una sucesión de CMTC, con po = (l/Z,l/Z), donde para

cada entero positivo n se t iene el generador infinitesimal An def inido

anteriormente.

Constrúyase la sucesión de PTUs {W(t,n)>:

t

O W(t,n) =S v(t,n)ds, W(0,n) = O para n=1,2 ,...

Entonces e’l PTU converge uniformemente a MB cas i seguramente, es de-

c i r:

P t l i m max IW(t,n) - w ( t ) l = o } = 1 n- O<t<T

INTERPRETACION DE LOS PTU.

Para cada entero positivo n, l a CMTC. v(t,n) es un proceso que t o - ma los va lo res -n y tn alternativamente, Inicialmente el proceso toma

ambos valores con igual probabilidad (l/Z), pero posteriormente va cambian

do de estados en los instantes a leator ios T1, T2, T3, ... De acuerdo con

el s iguiente teorema que ha s ido demostrado por Doob,J.L. (1953)], e l - tiempo de s a l i d a de cada estado tiene una distribución exponenci al con pa - rámetro n.

TEOREMA 2.

Sea A = . ) el generador inf in ites imal de l a CMTC v(t,n) y (ai ,J supongamos que

66.

para i = l , 2 , . p . ,N .

Dado v ( O ) , defini,mos el tiempo de sal ida como

T = in f {t > 0 I v ( t ) + i} Entonces :

P{T > t I v ( 0 ) = i l = 4 1 , i = l , 2 , 0 . . , N -x. t

En nuestro caso, tenemos sólo 2 estados (N=2), donde

x1 = x2 = n 2

as í que el tiempo de salida de cualquiera de los dos estados se distribuye según:

P{Ti - T i - l > t I v ( T i e l ) = 2 n } = e - n 2 t para i = 1 , 2 , ... El teorema 1 nos asegura 1 a aproximación de 1 as trayectorias del

PTU a una de las del MB, mientras que el teorema 2 nos dice cómo generar 1 os intervalos de tiempo en que el PTU se mantiene con velocidad consta; t e al ternante +n, para u n val o r dado de n .

En cuanto a 1 a rapidez de convergencia de los PTU. W hacia el MB W; [Vazquez Abad, F. (1984)], probó las hipótesis de independencia y normalidad de incrementos de MB simulado por PTU, utilizando la prueba j i cuadrada y encontró que tales hipótesis dejaban de rechazarse para valores de n ta les que

n > L-1 donde DT = longitud de intervalo de muestre0 de la t rayector ia .

67.

Debi'do a este resultado, podemos generar MB mediante PTU u t i l i zando un valor de n que sa t i s f aga l a desi'gualdad anterior.

3.5 Verificación posterior

Como se vió en el capítulo anterior, a cada trayectoria de MB simulado l e corresponderá una trayectoria solución de l a EDE ( l . l ) , l a cual tendrá barreras absorbentes en

x = o y X = l

De acuerdo con lo an te r ior , nos interesa que e l MB simulah ten - ga u n buen comportamiento en cuanto al alcance de barreras absorbmtes, para que no vaya a introducir u n sesgo en la trayectoria solución de (1.1).

Consideraydo un intervalo a,b , si el punto de partida del MB es x (con a < x < b ) , como se vió en el ejemplo 1.1.1, l a probabilidad de absorción de las trayectorias por 1 a barrera X = a es:

P(x,a) = (b-x)/(b-a)

y por 1 a barrera X = b es:

P(x ,b) = (x -a)/(b-a)

Entonces, para probar el buen comportamiento de las t rayector ias de M6 en cuanto a la absorción por las barreras X = O y X = 1 , bas- t a r á con probar que las proporciones de trayectorias que son absorbidas por ambas barreras, satisfacen (1 ).

En este caso hay que probar la hipótesis nula:

Ho : P(x , l ) = x

68.

contra la hipótesis al ternati.ya

Hi:P(x,l) + x

Como l a suma de absorciones por X = l , se distribuye binomialmente con parámetros N y x , entonces por l a aproximaci6n binomial a 1 a normal, tenemos que cuando N ti'ende a i n f i n i t o , l a proporción de absorciones - por X = l , tiene una distribución asintdtica normal , con media x y va - r i anza x(1-x)/N. De esta forma, para N grande , l a región de no rechazo de l a prueba de la hipótesis Ho anterior, con un nivel de significan - ci a a = O. 5 será [ LINF, LSUP) , donde :

LINF = X .. 1.964-

LSUP = x + 1.964'

El programa ESTO/ESlADISTICAS/Z que se presenta en e l APENDICE A , obtiene l a proporción de trayectorias que se absorben por cada barrera, a s í como los límites de la r e g i ó n de no rechazo, l o cual permite efectuar l a prueba para valores deseados de b-a y N.

69.

RESULTADOS DE LA PRUEBA, GENERACION MEDIANTF C4

romando en cuenta que, como ya hemos v i s to , para 1 as trayectorias de MB simuladas medi.ante CA se requiere que e l nümero de sa l tos de l a caminata por cada intervalo de muestreo , sea mayor que 40, se corri 6 el - programa ESlO/ESTADISTICAS/2 para N=100, a=O, b = l , y con nueve condi c io - nes in i c i a l e s , X = O . l , 0.2, . .. ,0.9.

Considerando 50 sa l to s de l a CA por cada intervalo de muestreo del MB, se encontró en todos los casos, que la hipótesis Ho no fue rechazada.

Esto se debe a que al simular trayectorias de MB mediante C A , se sat isface la hipótesis Ho, debido a que e s t a hipatesis también es sa- , i d ,"echa por las C A B +"..

GENERACION MEDIANTE PTU

Para probar la h ipótes i s Ho, mencionada anteri.ormente, se corr ió el programa ESTO/ESIADISTICAS/2, seleccionando l a generación de MB mediante PTU, dando On incremento de muestreo DT=0.001 y eligiendo un va - l o r igual a 75 para el parámetro velocidad, ya que según [Vazquez-Abad (OP. ci t 11 3

vel = [~W€I] = [ A m ] = [EZ?RJ'J = 75

Se encontró también en este caso, que para l a s 10 salidas elegidas de las t rayector ias (0.1,0.2, ..., 0.91, l a h ipótes i s Ho no fue rechazada.

De aquí se concluye que e l MB generado t isface las condiciones dadas en (1 1

. A continuación se presenta el l i s t ado de

en estas pruebas, tanto para la simulación por

por CA y por PTU, sa-

10s resultados obtenidos CA como por PTU.

70.

CONTEO DE SALLDAS POR BARRERAS SUP. E LNF, XO(1) : DIST. DE LA BARRERA SUP. AL PUNTO DE PARTIDA

PRI(1): PROP. DE SALIDAS POR LA BARRERA 1.NFERIOR PRS(1): PROP. DE SALLDAS POR LA BARRERA SUPERIOR

MB. SIMULADO POR: CA.

NUM. DE TRAYECTORIAS: 100 INTERVALO DE MUESTREO: 0.00190 VELOCIDAD: 50 NUM. DE BARRERAS SIMETRICAS EN CERO: 1 0

DISTANCIA ENTRE BARRERAS: 1.00 I XO(1) PRS(1) LINF(1) LSUPCI:]

1 0.10 2 0.20

3 0.30

4 0.40 5 0.50 6 0.60

7 0.70

8 0.80 9 0.90

O. 080 O. 180

O. 280

O. 440 O. 579 O. 650

O. 700

O. 810 O. 920

0.0412 O. 1216

o. 2102 O. 3040

O. 4020 O. 5040 O. 6102

O. 7216 O. 8412

O. 1588 O. 2784

O. 3898

O. 4960 O. 5980 O .6960

O. 7898 O. 8784 O. 9588

#ET = 7:27.7 PT = 2:54.5 I O = 0.1

CONTEO DE SALIDAS POR BARRERAS SUP. E LNF. XO(1) : DIST. DE LA BARRERA SUP. AL PUNTO DE PARTIDA PR I ( I ) : PROP. DE SAL IDAS POR LA BARRERA LNFERLOR PRS(1): PROP. DE SALI.DAS POR LA BARRERA SUPERIOR

MB. SIMULADO POR: PTU. NUM. DE TRAYECTORIAS: 100 INTERVALO DE MUESTREO: 0.00100 VELOCIDAD: 75

71.

NUM. DE BARRERAS SU4ETRI:CAS EN CERO: 10 DISTANCIA ENTRE BARRERAS: 1,OG. I XO( 11 PRS(I] LINF(I1 LSUPCI)

1 0.10 0.150 0.0412 0.1588 2 0.20 0.270 0.1216 0.2784 3 0.30 0.340 0,2102 0.3898 4 0.40 0.430 0.3040 0.4960 5 0.50 0.520 0.4020 0.5980 6 0.60 0.660 0.5040 0.6960 7 0.70 0.730 0.6102 0.7898 8 0.80 0.820 0,7216 0.8784 9 0.90 0.920 0.8412 0.9588

#ET = 2:23.2 PT = 1:Ol.g IO = 0.1

C A P I T U L O I v 1. Introducción.

En los capí tulos anter iores se ha presentado un método especifico que permite simular numéricamente la solución de ecuaciones diferenciales estocásticas haciendo uso de aproximaciones al proceso de Wiener basadas en caminatas a lea tor ias y procesos de transporte uniforme. En el presente capítulo se reseñará un enfoque al ternat ivo encaminado también a l a simu- lación de procesos de difusión pero que, a diferencia del anter ior , hace uso de cadenas de Markoff y de la discret ización de l a coordenada espacial , no l a temporal , para lograr este ob je t ivo . Para la presentación de este método nos basaremos fundamentalmente en I D i Masi & Runggaldier,(l981) I , referencia a l a que remitimos al lector interesado para un tratamiento mas detallado y extenso de lo que aquí se expone.

Considérese entonces l a ecuación diferencial estocástica en el sentido de I to

dxt = b(xt)dt + o(xt)dw x ( t ) = x

O O donde

b, o: [O, TI X Rn -> R n

son uniformemente Lipschitz y acotadas, w es un proceso Wiener estandar adaptado a l a f i l t r ac ión u sobre un espacio de probabilidad (Q,u,P).

La idea que subyace en el método que se presenta consiste en 1 a discretización de l a coordenada espacial de (1). Se obtiene así x , una Cadena de Markoff a Tiempo Continua (CMTC) con l a propiedad de que, al tender a cero l a norma de la d i scre t izac ión , se genera una sucesión de CMTC convergente al proceso de difusión solución de (1).

73.

Antes de proceder directamente a l a presentación de este enfoque se considera conveniente dar en ests sección algunos resultados previos que serán utilizados a lo largo del capítulo.

Sea entonces X := {Xt: t&TI una CMTC con T = [ O , l ) , E = { l , Z , . . . ,NI y = Z8, donde E es el espacio de estados con N E W y 8 es la algebra sobre l a que esta definida l a medida de probabilidad P. Las probabilidades de transición de l a CMTC X están dadas por:

P i j ( t ) = P ( X = j: X = i )

P . . ( O ) = "j 1J

que satisfacen como se sabe I Doob, 19531 :

( i ) min P . . ( t ) > O i 13 "N

( i i ) 1 Pi j ( t ) = 1 j=l

En consecuencia, l a matriz

P(t) := (Pij)

es una matriz estocástica con P(0) = I

( i i i ) P . . ( t + s ) = 1J

l a ecuación ( i i i ) , llamada notación matri c i al como

P ( t + S ) = P(s)P(t)

En consecuencia l a grupo es tocás t i co asociado

N

k=O 1 Pik(s)Pkj ( t ) t , s > o

i , j = 1 , ..., N

de Chapman Kolmogoroff puede reescribi rse en

familia P := { P ( t ) , t 2 0) constituye un semi- a la CMTC.

74.

El teorema siguiente asegura que el semigrupo correspondiente a una CMTC -en general a cualquier proceso de Markoff-, la caracteriza en el sentido explicitado en el Teorema, hecho que permitirá, posteriormen- t e , j u s t i f i c a r e l método de aproximación que aquí presentamos. Sea P el vector P(X = i ) , i = l , Z , ... ,N llamado la distribución inicial de X.

TEOREMA (Kurtz, 1981)

Sea X. un proceso de Markoff con espacio de estados E y distribu- ción inicial P. Sea P.:= P ( t ) el semigrupo correspondiente a X. .Entonces P. y P determinan las distribuciones de dimensión f i n i t a de X.

Además s i suponemos que I? es diferenciable en t = O , puede demostrarse que P. es diferenciable para todo t - > O y , en consecuencia si definimos A:= P(0)

b = AP

P = PA

P(0) = I

P(0) = I

y entonces

P(t) = exp (At)

La matriz A se l lana el generador infinitesimal del semigrupo, pues lo determina univocamente.

Sea ahora A:= (a . .). Se satisfacen entonces las siguientes condi- 1J

ciones:

min ai > O i = ly...yN j # i

N .. 1 ai j = 0 j= 1

75 .

Ahora bien, sabiendo ya que los semigrupos estocásticos determinan sus procesos de Markoff, es un hecho conocido I Hernandez, 19851 que s i construimos una sucesión de procesos de Markoff x que aproximen a o t r o , entonces

s i Ail -> A adecuadamente.

donde A y AN son los operadores infinitesimales respectivos de x y x . El sentido en que AN -> A puede precisarse de la s iguiente manera. Sean T y rN operadores lineales de D(A) en D(AN) los dominios respectivos de A y AN.

N

Entonces debe ocurr i r que

cuando N -> O0 con

Se entiende también, en el contexto de los procesos de Markoff, que l a condición de convergencia de los semigrupos respectivos significa que .

PN -> P s i

uniformemente para t l 0 ,T I y con + medible y acotada sobre E. Reescri- biendo (4) dicha 1 a condición de convergencia se convierte en

76.

para I) ~ C 1 5 1 1 dada por +(w) = $ ( w ( t ) ) IKurtz, 1981; Hernández, 19851.

En la siguiente sección se aplicarán los resultados anteriores para el establecimiento del método de aproximación.

2. Método de d i rcretización. -

E l problema que nos ocupará de aquí en adel ante es el de hallar generadores infini tesimales que sean "buenas aproximaciones" del generador infinitesimal del proceso de difusión solución de ( 1 ) y que garanti- cen l a convergencia en el sentido de ( 5 ) .

Con este fin en mente considérese que b y u son ta les que pasa cada condición in ic ia l x , ex i s te una solución única de (1) .

Sea E el espacio de estados definido previamente y tómese n = 1 puesto que se trabajará solamente en una dimensión por comodidad.

Sea h = T/N

Definase l a matriz A := como sigue h

A h = o i j

donde

77 .

es facil ver que

luego l a matriz

sat is face (3) de tal forma que constituye el generador infinitesimal del semigrupo de una CMTC {xt: O < t < TI con espacio de estados E = {1 ,2 , . .. , N ) que sat is face

h "

P(x ( O ) = x 1 = 1 h O

Dado que AN cumple con ( 6 ) en particular ocurre que

Ahora bién, hasta aquí hemos obtenido en general una posible aproximación para l a solución de (1) que, sin embargo, no toma en cuenta l a necesidad práctica de traba jar en regiones acotadas. Esto obliga a tomar en cuenta los tiempos de primera salida de l a región y el proceso parado asociado que definimos de 1 a forma siguiente:

Sea CCR abierto y acotado y sea Ch:= C n E. Entonces definimos el instante de primera salida de C mediante

78.

T x t s C

inf{t: ts para x.

s i t s l O , T I í T x t ~ C h t E

l 0 , T I xttC1 s i no Th = [ i n f ( t : ; E / O , T

para x .

Definimos los correspondientes procesos parados mediante

A

Xt - Xt A T .. -

Habiendonos restringido a una región acotada de R debemos especi- f i c a r l a s condiciones que debe s a t i s f a c e r l a CMTC en la f rontera de l a región, definida por

Es claro también que el espacio de estados de 2 , es C U 3C . h h

En lo que sigue consideraremos e l caso en el que tanto el estado inic ia l como el final de 1 a CMTC aproximada son absorbentes

Debe decirse, sin embargo, que esta elección es arbitraria y pu- ramente i lustrativa. Dependiendo de l a situación específica a modelar, l as condiciones a la frontera cambiarán. En particular en esta monografía se analizará el caso de un proceso de difusión proveniente de la genética de poblaciones en el que se consideran estados absorbentes tanto el estado final como el in ic ia l .

3. Ccnvergencia del método.

Para lo que sigue nos basaremos fundamentalmente en el desarrollo que hacen I Di Masi & Runggaldier (1981) I bajo las condiciones (8) y deja-

79.

remos para la sección última de esta monografía l a discusión de los resultados computacionales obtenidos.

Se da a ccntinuación una demostración de l a convergencia de una sucesión de CMTC a l a solución de ( 1 ) . Para tal efecto denótese por D l 0 , T I al espacio de funciones que toman valores en Rn sobre l0 ,TI , continuas por 1 a derecha, continuas por l a izquierda en T y con 1 ímites iz- quierdos equipado con l a métrica de Skorojod con el objeto de que D n l O , T I resulte separable y completo.

n

Debido a que dar un proceso estocástico sobre un espacio probabi- 1 i s t i co , equivale a proporcionar una p-robabi 1 idad sobre el espacio de 1 as trayectorias del mismo, una sucesión de CMTC genera una sucesión {PN} de probabilidades sobre el espacio de 1 as trayectorias de 1 a cadena , d i - gamos (Q,a); en consecuencia, s i deseamos que PN -> P cuando N -> 00

debe ocurrir que

1 im /$ dPN = /$dP N- Q S1

para todo $eClQl

I Di Masi & Runggal d i e r , op. c i t . 1 demuestran que las cadenas generadas según el algoritmo dado en l a sección 2 convergen debilmente a la solu- ción de (1) . Este tipo de convergencia basta puesto que estamos interesados solamente a l a convergencia de funcionales de las trayectorias del proceso estocástico considerado.

Considérese ahora una familia {x 1 de CMTC que aproximan l a d i f u - h

sión original x. Para cada h > O , definamos x por h

xh = x + Bh( t ) + l h ( t ) O

con

80.

Bh(t ) = Jt O b(x!)ds

Los autores citados demuestran que las sucesiones B ( t ) y 1 ( t ) poseen subsucesiones convergentes respectivas, basándose en el hecho de que la existencia de estas subsucesiones es equivalente a la t i rantez (en inglés tightness) de l a sucesión original . Mas concretamente s e t i e - ne que

h h

se sat is face para (9b) y (9c) para c , g , K PO s i t ivos y co Insta mtes. Esta condición garantiza la t irantez, luego, existen subsucesiones convergentes IBil l ingsley, op. c i t . I .

El problema que resta es el de l a demostración de la ex i s tenc ia de un Gnico l imite para las subsucesiones que resyl ta ser x . , l a solución de (1).

La idea consiste en demostrar que e l proceso l ímite está dado por

donde w^ es un proceso estandar de Wiener adecuado (no necesariamente el original en (1) . Nótese que se ha relajado un tanto el concepto de solu- ción de (1) ) .

En realidad l o que se t ra ta es de ha1 l a r e l adecuado para que x converja a x. y , por lo tanto , debe demostrarse que B ( t ) -> B ( t ) . h h

81.

No nos detendremos a comentar los detalles técnicos de esta par te de l a demostración ya que, por su naturaleza, pensamos que el lector interesado debe remit i rse a l ar t ículo or iginal .

Se tiene entonces demostrada 1 a corivergenci a débil de 1 a CMTC a l a solución de (1). Para el proceso parado ct ocurre exactamente l o mismo, es decir , xt convcrgc a xt. Debe recordarse que e s t e proceso parado es el que se usará para l a implementaci6n práctica del nétodo propuesto , cuya mecánica es como sigue:

Ah A

i ) Se obtiene la matriz A := ( X . .) de orden n x n que se i n - h h h h 1 SJ

crernenta al disminuir el psso de discretización h . Sea ahora

el instante en que se dá el primer s a l t o en l a CMTC x.

Por el teorema 2 del capítulo 3 tenemos que

don de

h h hi:= xi ,i+l + Xi

Más aún, los instantes en que se dan los saltos sucesivos forman una sucesión monótona creciente

donde

i = 1, 2, 3, ...

82.

Así pues, basta saber generar variables exponenciales con pará- metros X , . . . ¶ X,,, para que podamos generar las CMTC respectivas ¶ cuando menos en lo que se ref iere a los instantes en que se cambia de estado. En cuanto a l a magnitud de los sal tos, hay que muestrear 1 a distribución dada por l a i-ésima hi lera de P ( t ) para obtener XT cuando X = i ; s i X = j , entonces se mantiene constante l a t r a y e c d r i a en j hasta el ins- t a k e T_ = T, + o , donde

1

O

T

L I

P(a > t : X = j) O

En ese instante

= exp( -x .t) J

se cambia al nuevo estado digamos R, etc.

Así pues , hay que conocer l a función matricial P(t) para lo cual prescribimos

i i ) P h ( t ) = Ah(t)Ph(t )

P ( O ) = I h

y P (t) es l a matriz de t r a n s i c i h de una Cadena de Markoff que aproxima x. cuando h es suficientemente pequeña.

h

El problema práctico que se presenta para l a implementación del método radica en l a solución numérica de las ecuaciones del paso i i ) . Da- do que estamos discretizando un operador de difusión con espectro continuo e infinito, los valores propios que se obtienen para i i ) difieren en valor absoluto grandemente, lo que l o convierte en un sistema de ecuaciones diferenciales r í g i d o .

Aunada a esta dificultad se encuentra otra que es I a relativa a l a determinación del tamaño adecuado de h , e l paso de discretización que asegure l a convergencia de x a x. h

C A P I T U L O V

1. Factores reguladores .de 1 a variabi 1 i d a d gerléti ca.

En es te capí tulo se proporciona al lector la justif icación bioló- gica del modelo par t icular del que se ha hecho uso en los capítulos an- t e r io re s para i lus t ra r l as t écn icas de simulación analizadas para EDE. Se comparan también los resultados numéricos obtenidos aquí con los resultados teóricos publicados, cuando esto es posible. Cabe, sin embargo, hacer l a aclaración de aue el objetivo de lo que sigue a continuación no es e l de presentar como nuevos resultados ya conocidos desde hace tiempo sino, por e1 contrario, unicamente el de i l u s t r a r l a bondad de las técnicas de simulación para EDE que se han mostrado en capítulos anteriores. Hecha es- t a aclaración procedemos al planteamiento de nuestro modelo.

PL;ede decirse que existen tres grupos de factores que cambian 1 as frecuencias génicas en las poblaciones de especies biológicas ICrow & Kimura, 1970 I .

E l primero de estos comprende l a s llamadas presiones evolutivas sistemáticas en las que están incluidas las tasas constantes de mutación y migración y la selección natural de intensidad constante.

E l segundo de los grupos consiste de los factores que producen fluctuaciones aleatorias y que pueden dividirse a su vez, en dos t ipos

a ) el llamado muestre0 aleator io de gametos en poblaciones f i - ni tas o deriva Sgnica.

b ) las f luctuaciones aleatorias en las presiones sistemáticas, por ejemplo, en el coeficiente selectivo.

84.

Finalmente tenemos a l grupo en el que se encuentran incluidos los rearregl os cromosómicos , dupl i caciones o def i ciencias de nucl eóti dos , po- l ip lodía , e tc .

E l modelo que estudiaremos en esta sección se centra en el segundo grupo de factores, particularmente los incluidos en el inciso a).

E l problema que el modelo puede ayudar a resolver es el de prede- c i r l a cantidad de polimorfismo y heterocidad que es mantenida en una po- blación de tamaño f i n i t o debida a l a acción conjunta de la selección natural y l a deriva génica.

Antes de proceder directamente al estudio del modelo, nos detendremos a comentar y proporcionar los conceptos y definiciones básicas de las que haremos uso. Los lectores para los cuales n G sea necesario esta introducción, pueden pasar directamente a l a siguiente sección que es i n - dependiente de ésta .

Se entenderá población er: este capítulo, como un grupo de i n d i v i - duos en los que es posible el entrecruzamiento y que coexisten en tiempo y espacio. Esta definición es sufScientemente operativa para los usos que aquí l e daremos y no pretende comprender en el 1 a todas las posibles acep- ciones del término.

Así mismo, se supondrá que existen dos sexos y que l a población está consti tuida de individuos sexualmente maduros que se aparean aleatoriamente i . e . : es un grupo en el cual l a probabilidad de apareamiento en- t r e individuos de genotipos particulares o fenotipos es igual al producto de sus frecuencias individuales en l a población. Este tipo de poblaciones se conoce también con el nombre de panmícticas.

Ahora bien, ya que nuestro interés es el de descr ibir la población en términos genéticos, debe examinarse la naturaleza genética de los caracteres fonotípicos bajo las restricciones determinadas por l a definición

85.

de población que se ha tomado aquí.

Considérese ent,onces e l caso de ccn locus autosómico con dos alelos A y a. El número de individuos de cualesquiera de los tres genotipos posibles AA, Aa y aa se corresponde con el nümero de individuos de las clases fenotípicas posibles cuando no exis te dominancia por parte de ninguno de los alelos. S i denotamos por n , n los individuos con genotipos A A , Aa y aa respectivamente ta les que

y n 2 2 1 1 1 2

n + n + n = N I 1 1 2 2 2

donde N es el número total de individuos en l a población, entonces definimos las frecuencias genotípicas por

D = n /N 3 1

H = n /N 1 2

R = n /N 2 2

donde, D , H y R representan las frecuencias relativas de AA, Aa y aa res- pecti vamen te.

Las frecuencias genotipicas indican entonces l a forma en que l a variación genética de l a población esta distribuida entre los diferentes genotipos.

Definanse ahora las frecuencias génicas como las razones

p = (2D + H)/2N

q = (2R + H)/2N

en donde p es 1 a fracción de alelos A en e l acervo genéti co y q es 1 a

86

fracción de alelos a.. Debe tenerse en cuenta que se está suponiendo que el número t o t a l de individuos diploides N se forma debido a l a unión de 2N gametos provenientes de l a generación anterior cuando el número to ta l de individuos no cambia.

Las frecuencias géni cas nos proporcionan información sobre l a va- riación genética dentro del acervo genético t o t a l de la población.

Ahora bien, en ausencia de factores que cambien las frecuencias génicas como la Selección, la mutación, la deriva génica o l a migración y bajo l a acción de apareamiento aleatorio, las frecuencias genotípicas de un locus simple después de una generación pueden representarse median- t e una distribución binomial relacionada con las frecuencias génicas.

E l enunciado anterior se conoce con el nombre del principio de Hardy-Wienberg y predice la existencia de un punto de equilibrio para las frecuencias genotípicas cuando las frecuencias génicas no cambian:

D = p2

H = 2pq

R = q2

con p2 + 2pq + q 2 = 1.

S i n embargo, en l a s poblaciones naturales dificilmente existen casos en los que no se presenten factores que al teren el valor de 1 as frecuencias génicas. En consecuencia, estamos intgresados en estudiar las ca- racter ís t icas par t iculares que producen sobre dichas frecuencias los factores rnerxionados.

Este trabajo se centrará fundamentalmente sobre dos de e l l o s : l a selección natural y la deriva génica.

87.

Para la caracterización del primero de los conceptos mencionados y con mucho e l mas importante del pensamiento evolutivo contemporáneo re- curramos a Darwin:

"...Las variaciones por l igeras que sean, y cualquiera que sea l a causa que las produzca, si son en algún grado provechosas a los i n d i v i - duos de una especie en sus relaciones.. .con otros seres orgánicos y en sus condiciones f í s i cas de vida,tenderán a l a conservación de estos individuos y serán en general heredadas por l a descendencia, l a cual a su vez también tendrá asi mayor probabilidad de sobrevivir. De lcs muchos individuos de una especie cualquiera que nacen periodicamente, solo un número pequeño puede sobrevivir. Este principio por el cual toda ltgera variación, si es ú t i l , se conserva, l o he denominado yo con el término de selección natural . . . ' l . "(cf . El Origen de las especies)?

Esquemáticamente y con terminología modern? podemos representar l a selección natural de la s igu ien te forma (Hedrick, 1984):

Imagínese una población con una cantidad dada de variación genéti- ca entre sus individuos. La variación genética en combinación con componentes ambientales, determina la var iación femtípica para e l conjunto de caracteres expresados por un organismo particular. La cantidad y naturaleza de 1 a variación fenotípica a su vez da 1 ugar a variaciones en sobrevi- vencia, fecundidad, habilidad de apareamiento, e tc . , fac tores que f i n a l - mente determinan l a transmisión hereditaria de alelos a generaciones posteriores.

VARlAClON GENOTIPICA AMBIENTAL FENOTIPICA EN EFICACIA

* ""~""""""""""~""""-""~""""""""""""""""~

t

Darwin C. (1859) E l Origen de las Especies. CONACYT, México.

88.

Una vez caracterizado el concepto de selección natural , se proce- derá a la presentación de l a terminología básica con l a que será posible después construir modelos que nos den 1 a posibi 1 idad de predecir trayectorias evolutivas de l a s poblaciones.

Defínase en consecuencia la e f icac ia re la t iva como la habilidad relat iva de diferentes genotipos para transmitirse a través de sus aielos a generaciones posteriores.

Considérese ahora la si tuación en l a que los genotipos AA, Aa y aa poseen eficacias relativas respectivas w , w , w y supongamos que q es la f recuencia a le l ica de a y 1-q l a de A .

11 1 2 2 2

La ef icacia re la t iva media 6 de l a población se define como l a suma de las contribuciones relativas de los diferentes genotipos. Antes de que actue la selección y por el principio de Hardy-Wiemberg, 1 a efica- c i a media está dada por I Hedrick31984j :

Obsérvese que 1 a contribución re1 ativa medi a de cada genotipo a l a siguiente generación se determina por el producto de la ef icacia re la- tiva respectiva y la frecuencia previa a l a acción de la selección del genotipo particular.

Si llamamos q a 1 a frecuencia de a después de l a acción de la se- 1

lección entonces

q l = (2pqw 1 2 /a/z + q 2 w 2 2 /i = (pqWl2 + q 2 w 2 2 )/i

trow & Kimura, 19701.

El valor del cambio en l a frecuencia génica l a definimos como

89.

1 uego

Aq = ( pqw12 + q2 w - qG)/d; 2 2

por lo tanto y del hecho que q = 1 - p obtenemos:

Observemos ahora que

s i $ es diferenciable. Luego

di/dq = 2[q(w - w ) - p(w 22 1 2 11 - WI2)1

y en consecuencia

Se tiene entonces que e l cambio en frecuencia a lé1 ica es función

de las f recuencias a lé l icas, ef icac ia media y el cambio en e f i cac i a media con respecto a l a f recuenc ia a lé l i ca .

Supondremos ahora que l a e f i c a c i a del heterocigoto es 1 y modifi-

cando adecuadamente l a s e f i c a c i a s de l o s homocigotos se estudiará el efec- t o que tiene wij sobre el cambio genét ico. S i se supone que w pueden ser d iferentes, existen cuatro t ipos de selección posible:

11 y "12

Selección por A

Selección por a

Ventaja del heterocigoto

90.

Desventaja del heterocigoto.

En par t icular nos centraremos sobre el caso de selección progre- s iva , en l a cual s e da la selección para alelos que son ventajosos para un ambiente dado.

En es ta s i tuación estamos suponiendo que l a población esta sujeta a cambios ambientales que pmdccen el reemplazamiento de alelos que den una mejor adaptación 2 través de la selección natural . Debe hacerse h i n - capié en que no estamos postul ando que nuevos alelos mutantes "surgen" a consecuencia de l o s c a h i o s ambientales ( i .e.: heredabilidad de los caracteres adquiridos), sino que alelos ya existentes en el acervo genético de 1 a población se ven favorecidos por e l medio y , en consecuencia, se selec- cionan en el 1 ocus al que pertenece.

Denótese por h el nivel de dominancia. La selección progresiva es función de h y por consiguiente, estamos interesados en estudiar e l cambio en l a frecuencia de A a p a r t i r de una frecuencia inicial . Tenemos entonces la siguiente si tuación a modelar:

G e n o t i p o

Adaptabilidad AA Aa aa W W W

l+s l+hs 1

con O " < S < 1 11 amado el coeficiente de selección. Luego susti tuyendo los valores del cuadro anter ior en (1) tenemos

donde d; = 1 + 2 h s x ( l - x ) + sx.

Para el caso S = 0.1 y x = . O 1 tenemos que A se f i j a en aproxima-

91.

damente 100 generaciones cuando no existe dominancia para a , es decir h = 1/2 (ver gráfica 1)

A I "

0.75 --

U t

050 --

o25 --

I I I I I I I

do 2bo 3bo &o 500 600 700

Gráfica 1: Tiempo de f i jac ión para A sin dominancia y para S = 0.1 . ..

x = 0.01 (Según Hedrick,l984). O

_. .

Terminamos aqui los conceptos que se han considerado necesarios para el entendimiento del modelo general de la siguiente sección y se ini- c ia ahora una breve discusión sobre l a deriva génica y su acción sobre l a variación genética.

La deriva génica puede definirse como el cambio que sufre la es- tructura del acervo génico de una población debido a efectos aleatorios que se producen durante 1 a transmisión del acervo génico de generación en generación. Este cambio es errático, de naturaleza también aleatoris. y

aunque débil comparado con el efecto de l a selección natural, es no obstante, relevante cuando el tamaño de l a población es pequeño [Roughgarden , 19 791 .

Es precisamente l a deriva génica el factor que introduce elemen-

92.

tos estocást icos en l a model ación de la va r iab i l idad genét ica de l a s po-

bl aciones.

Para estudiar el efecto de la der i va gén ica en poblaciones pe-

quefias se cuenta con el modelo de Wrigth-Fisher mediante el cual se puede ca l cu l a r l a p rob sb i l i dad de que 1 a frecuencia génica tome a lgún va lor

preestablecido. El modelo es una cadena de Msrkoff con probabil idades de

t rans i c i ón dadas po r

y con d i s t . r i buc i ón i n i c i a l dada.

Estas probabi l idades de t rans ic ión representan la probabi l idad de l a e x i s t enc i a de i a le l o s a en l a generación t+l dado que e x i s t í a n j de e l l o s en l a generación t.

En l a g r á f i c a (2) se presenta una simulación, haciendo uso de (2),

de l a d i s t r i b u c i ón de las f recuencias a lé l icas para una población de tama-

ño N = 20.

1 x) GENERACIONES

Gráfica 2: Simulación del efecto de la der iva génica sobre una población de tamaño N = 20, x. = .5 pa r a l a lo. ,50. y 20. generación.

93.

Puede observarse que el efecto final de la deriva génica es el de convertir en homocigota para alguno de los dos alelGs ;1 l a población. E’! proceso es sumamente ráp ido y eventualmente l a m i tad de l a población se f i j a r á para a y l a o t r a para A sin existencia de heterocigotos.

Por otro lado, puede calcularse también la dis t r ibución de las frecuencias de f i jación para e l a le lo a por ejemplc, haciendo uso de l a matriz de transición de Wright-Fisher. En l a f i g u r a (3) se muestran estas distribuciones (Hedrick, 1984) para diversos valores de la frecuencia i n i ci al .

Puede observarse que conforme el valor de 1 - x decrece, la moda de l a distribución se desplaza hacia la derecha. Este fenómeno es ref le jo de la canti dad de cambio total requerido en 1 as frecuencias alé1 i cas para obtener la fijación final del a le lo .

7%

a04 -- q = O 0.8

O 50

GENERACIONES

200

Gráfica 3: Distribuciones de las poblaciones fijadas para a para las frecuencias iniciales x = .8, .5 y .2 y N = 20.

94.

2 . Modelo de d i f u s i ó n para seimción n a t u r a l y deriva génica.

El problema que nos ocupa en esta sección es el de predecir 1 a cantidad de polimorfisno y heterocigocidad que se puede mantener en una población de tamaño pequeño cuando es ta sometida a selección progresiva y deriva génica.

Para poder hacerlo aproximaremos el proceso de cambio en las f recuencias génicas como un proceso estocástico, en par t icular como u n proceso de difusión [Crow 6 Kimura, 1970; Kimura & Otha, 19711.

Este enfoque de modelación nos permi tira’ después , apl icar l as herramientas desarrolladas en es ta monografía para e l anál is is del modelo resultante.

Considérese entonces la si tuación en que un pa r de alelos A y a se segregan en u n a población de tamaño N con frecuencias respectivas x y 1 - x. Debe suponerse también que A es u n a le lo mutante de a que es favo- recido cuando las condiciones ambientales lo propician y que t iene una l igera ventaja selectiva S sobre a en esas situaciones. Se supone también que l a población es panmíctica. Estamos entonces ante un fenómeno de se- lección natural progresiva.

E l objeto que se persigue en lo que resta de l a sección es l a de encontrar una expresión para l a probabilidad de f i jación de A en l a t-ésima generación dado que su frecuencia inicial es x en t = O .

O

Enfocaremos , como ya se ha dicho , 1 a model acjór: de es t e fenómeno suponiendo que e l proceso de cambio en la frecuencia génica de A es un proceso de difusión. La justificación rigurosa de esta hipótesis puede encontrarse en [Kimura, 1970; Crow & Kimura, 1970 y Kimura & Otha , 19711, por ejemplo.

Recurriendo a l a fórmula (1) tomando h = 1/2 y suponiendo S muy

95.

considerar como valores a la frontera

x( .02) = o ~ ( - 9 8 ) = 1

Ambas ccndiciones (4a . ) y (4b.) son de l a forma

x(a) = O x(b) = 1 4

y significan que l a probabilidad de que el alelo A s e f i j e cuando x = a es nul a y que l a misma probabi 1 idad es 1 cuando x = b.

O

O

Con las condiciones a la frontera ( 4 a . ) l a solución de (2.4.4.) es

O

donde u representa l a probabilidad última de f i jación dado que l a frecuencia inicial es x .

O

Para e l caso resuelto en las simulaciones del capítulo 2 debe re- solverse (2.4.4.) con las condiciones (4b) obteniendose que

u ( x o ) =

es 1 a expresión

exp(-4Ns( .02)) - exp(-4Nsx ) exp(-4Ns(.02)) -exp(-4Ns(.98))

anal í t ica que se aproxima con l a simulación a través de los siguientes tiempos de paro:

en donde T e s e l estimador de l a probabilidad de que e l a l e l o s e f i j e en l a población y T e s e l estimador correspondiente de l a probabilidad de desaparición.

1

2

96.

pequeña y haciendo las sustituciones pertinentes para x la frecuencia de A. obtenemos que

luego el coeficiente de deriva es

b(xt) = (s/Z)[x( 1 - x)] [Crow & Kimura, op. cit.]

E n cuanto al coeficiente de difusión [Crow & Kimura, op. cit.] demuestran que está dado por l a vari anza de la distribución binomial de l as f recuan ci as al é1 i cas

a(Xt) = ~ ( l - x)/2N

con estas funciones así definidas la EDE resultante es

dxt = sx(1 - x)dt + (x(1 - x)/2N)’/‘dw

x(0) = x O

donde hemos reetiquetado S : = s/2.

Resumiendo, tenemos que ( 3 ) es el modelo de selección progresiva de u n alelo mutante A con ventaja selectiva S en una población panmíctica y bajo el efecto de ? a deriva génica.

A esta ecuación esta asociada l a ecuación (2.4.4. ) con las condiciones a l a frontera

x ( 0 ) = o x(1) = 1 (4a. )

Sin embargo, l a ecuación ( 3 ) presenta problemas de rigidez en l a frontera al resolverla numéricamente, l o que conlleva a l a necesidad de

97.

Las correspondientes estimaciones para ( 5 ) estan dadas por los tiempos de paro

No estudiaremos aqui en l a presentación de los resultados analí t i- cos conocidos y publicados para el modelo con el que trabajamos ahora. Al lector interesado en p r o f u n d i z a r en éstos y otros aspectos de la genética de poblaciones , l e recomendamos I Crow & Kimura, 19701 amp1 iamente.

En la sección siguiente daremos la interpretación de los resultados de las simulaciones presentadas en el capitulo 2 y compararemos en a l - gunos casos , estos resultados con los teóricos cuyas fórmulas hemos presentado en esta sección.

3. Disaasión de resul tadar.

Iniciaremos esta sección especificando los valores de ICs paráme- tras de (3) con los cuales se trabajará.

Considerando que (3) model a el efecto de 1 a deriva génica y 1 a se- lección natural tiene sobre la heterocigosidad y polimorfismo de un locus en una población de tamaño pequeño , se ha tomado N = 10 con el objeto de apreciar mas claramente el efecto del coeficiente de deriva de (3) sobre l a dinámica evolutiva de l a población. Cosa s imilar sucede con el coeficiente selectivo s. Se ha tomado aquí como S = 0.1 como un valor repre- sentativo de ventajas selectivas l igeras del a le lo mutante.

En cuanto a l a condición in i c i a l , s e han tomado diversos valores de x que son

O

x = .1, .2, .3, .4, . 5 , .7 , .8 O

98.

En el cuadro (1) se presentan los resultados correspondientes a las probabilidades de f i jación según las fórmulas ( 5a . ) y (5b.) y los tiempos de paro obtenidos p o r simulación de ( 3 ) para cada condición i n i - c ia l x . Se demuestra también el error obtenido entre el valor real dado p o r (5b.) y la simulación.

CUADRO 1: Cuadro comparativo de las probabilidades de fijación obtenidas por simulación y 1 as fórmulss ( 5 ) .

X 5a. 5 b . Sim. Error

.1 .3358

.2 .5609

.3 ,7122

.4 .813

.5 .8808

.6 .9262

- 7 .9562

- 8 .9771

.2 799

.5245

.6885

.7984

.872 1

.9768

.9546

.9 768

.28

.42

.72

.78

.82

.91

.95

1

. O00 1 -1045

- .O315

. O 184

.O5 10

.O115

.O046

- .O232

En l a g ráf ica ( 4 ) puede apreciarse también de los resultados simu- 1 ados .

1.00 " c

o. 75 "

0.50 "

0.25 "

I I I I 1 I I

0.1 0.2 0.3 0.4 0.5 0.6 0.7 08 0.9 1.0

Gráfica 4: Se muestran los resultados obtenidos con l a fórmula (5b.) y l a simulación de (3) .

99 .

Observamos entonces que la probabilidad de fi jación de A aumenta conforme la frecuencia inicial es mas elevada. Este es u n resul tad0 espe- rado ya que estzmos considerando que A t iene una ventaja selectiva sobre a , aspecto que favorece su f i jación. Puede compararse este resultado con la gráfica 2 en donde consideramos el efecto aislado de la selección natu ra l . E l alelo A termina fijandose en alrededor de 100 generaciones. La seiección actúa favoreciendo e¡ alelo que l e concede a l a población mayor posibilidad de sobrevivir en el ambiente cambiante.

En l a sucesión de gráficas siguientes se presenta la distribución de las Probabilidades de f i jación para A cuyos valores se han obtenido promediando los correspondientes obtenidos en el capítulo 2 a través de Caminatas Aleatorias y P.T.U.

Podemos observar el mismo fenómeno: mientras mas elevada sea la frecuencia inicial de A en la población, mayor es la probabilidad de su f i j a c i i n .

x = 0.1

40.5

X = 0.4

57

0.0 02 0.4 0.6 08 QO 0.2 Q4 Q6 0.8 0.0 1

X = 0.80

10.5 -

0.2 0.4 Q6 Q8

Gráfica 6: Distribución de las probabilidades de fijación para diversos valores de x (Promedio de valores obtenidos por Caminatas Aleatorias y P.T.U.).

O

100.

F i n a l m e n t e p r e s e n t a m o s e n l a g r á f i c a ( 7 ) ¡as d i s t r i b u c i o n e s de l a

f r e c u e n c i a g é n i c a de A p a r a d i f e r e n t e s g e n e r a c i o n e s c o n s i d e r a n d o u n a f r e -

c u e n c i a i n i c i a l x = .3. En es tas cond ic iones es tamos supon iendo que tan

s o l o e l 33% de l a p o b l a c i ó n c u e n t a , y a s e a h o m o c i g ó t i c a o h e t e r o c i g ó t i c a -

mente, con e l a l e l o A en su g e n o t i p o .

O

t = 5

175 125 135

O.@ 02 0.4 0.6 0.8 1.0

t = 20

235 f

t = IO t = 15

t - 3 5 55.5 t =50

28

'i 45 2.

I I 0.5

ao a2 0.4 0.6 08 1.0 0.0 02 0.4 0.6 a8

29

I 0.5 0.5 0.5

68

1.0 0.0 0.2 0.4 0.6 0.8 1.0

G r á f i c a 7: D i s t r i b u c i ó n de l a f r e c u e n c i a g é n i c a de A en 100 p o b l a c i o n e s

p a r a d i f e r e n t e s t i e m p o s d o n d e x = .3. O

101.

Es a q u í donde puede apreciarse el efecto combinado de la deriva y la selección. Si comparamos estas gráficas con las gráficas correspondientes a los efectos de la selección y deriva aislados sobre 1 as probabilidades de f i jación. Observamos que la deriva tiende a volver homocigo- tas a las poblaciones para u n o u o t r o alelo, es decir elimina la variabilidad genética. Sin embargo, cuando la selección actúa también , las cosas cambian. Existen siempre poblaciones que contarán con individuos polimór- ficos con genotipo Aa y será una mayor proporción de poblaciones l a que se f i j a rá para A comparada con l a proporción que s e f i j e para a .

La selección natural y la deriva génica son entonces fuerzas evolutivas de naturaleza opuesta. La segunda de el las t iende a homogenizar e l acervo génico de 1 as poblaciones y l a primera l e da variabilidad que en términos evolutivos significa mayor adaptabilidad y , en consecuencia, mayor capacidad de sobrevi vi r.

Debe indicarse, no obstante, que l a deriva génica solo puede se r s ignif icat iva p a r a poblaciones de tamaño pequeño ya que , según puede apreciarse de (3) , al aumentar el tamaño de l a población N el coeficien- t e de difusión tiende a cero y, en consecuencia, para poblaciones grandes, solo contará el efecto de la selección natural.

Un comentario resta por hacer. Las técnicas de simulación son en general herramientas versátiles y con una potencialidad de manejo para los científicos no matemáticos, muy grande debido a su r e l a t jva f ac i l i - dad de implementación y al uso de estimadores estadísticos con los que este t ipo de cient í f icos es ta famil iar izado, además de contar con paque- tes especializados para este t i p o de anál is is . Esperamos que el ejemplo que aquí se presenta motive su uso en mayor escala por parte de e s t e t i - po de cient í f icos .

A P E N D I C E A

A q u í se presentan los programas uti l izados en el capítulo 3 en l a generación y pruebas de M B . Aunque en el cuerpo de los programas aparece una documentación de los mismos , haremos una descripción complementaria y presentzremos alguna información auxi l ia r que f a c i l i t e su uso.

SUBRUTINAS.

CAMINA (SEM , NINC , INCR, DT, DCAM)

Esta subrutina simula incrementos de NB mediante CA; estos incrementos se utilizan posteriormente en otros programas. Para cada incremen- t o de MB se simulan [ DT/DCAM) sal tos de una CA. De acuerdo a las pruebas efectuadas, esta subrutina simula de manera aceptable incrementos de trayectorias de M B , cuando LIT y DCAM son ta les que:

[DT/DCAM] > 40

PTU(SEM, NINC, INCR, DT, N )

Esta subrutina genera incrementos de MB mediante PTU, los cuales se uti l izan después en otros programas. En la f igura 1 se tiene el diagra- ma de flujo uti l izado en la construcción de la subrutina, el cual fué tomado de [Vázquez-Abad, F. (op . ci t)] .

El parámetro N (velocidad) debe sa t i s f ace r , como ya hemos v is to , l a desigualdad:

103.

s i

FIGURA 1. Algoritmo para generar Movimiento Brownian0 mediante Procesos de Transporte Uniforme.

104.

L ESTO/MONO

# F I L E ( M E E O ) ESTO/PlONO ON PACK

10000

10100 10200

10300 SUBROUTINE PTU(SEM,NINC, INCR,DT,N)

10400

10500 C* SEM = S E P I I L L A

10600 C * N I N C = NUM. DE INCREMENTOS DEL PROCESO A GENERAR

10700 C * I N C R = ARREGLO UNIDIMENSIONAL y DONDE SE GUARDAN LOS INCREMENTOS

10800 C* GENERADOS

10900 C* DT = INCREMENTO DE TIEMPO DEL PROCESO GENERADO

11000 C* N = I N D I C E D E L T E R M I N O U T I L I Z A D O DE LA SUCESION DE PROCESOS (PTU)

11 100 11200 R E A L I N C R ( N I N C )

11300

11400 C * I N I C I A L I Z A C I O N

11500 J = O

11600 SANT = O

1 1700 s = o 11800 T 1 = O

11900 T2 = O

12000 V E L F L O A T ( N )

12100 IF(RANDOM(SEM) .LT. O. 5 ) V E L = -VEL

12200

12300 C* I N I C I A N I T E R A C I O N E S

12400 DO 100 I I N C = 1 ,NINC

12 500 J = J+l

12600 T I = I I N C * D T

12 700 I F ( TAU . GT. DT) GOT0 80

12800 60 S = S + VEL*TAU

12900 T1 = T 1 + T A U

13000 u = RANDOM(SEM)

13100 TAU = -ALOG(U)/(VEL*VEL) 13200 T2 = T2 + T A U

105.

13300 VEL = - V E L

13400 I F ( T 2 .LT. T I ) GOTO 60 13500 T = T I - T 1

13600 GOTO 90 13700 80 T = DT

13800 90 S = S + V E L * T

13900 TAU = TAU - T

14000 T 1 = T 1 + T

14100 C* SE GUARDAN INCREMENTCS DEL PROCESO

14200 I N C R ( J ) = S - SANT

14300 SANT = S

14400 100 CONTINUE

14500

14600 RETURN

14700

14800 END

14900

150CO

15100

15200

15300 SUBROUTINE CAMINA(SEPI,NINC,INCR,DT,DCAM) 15400

15500 C* SEM = SEMILLA , DEBE ESTAR EN ( O ,1) 15600 C * N I N C = NUM. DE INCREMENTOS DEL PROCESO A GENERAR

15700 C* I N C R = ARREGLO UNIDIMENSIONAL, DONDE SE GUARDAN LOS INCREMENTOS

15800 C * GEN E RADOS

15900 C* DT = INCREMENTO DE TIENPO DEL PROCESO GENERADO

16000 C* DCAM = INCREMENTO DE TIEMPO PARA EL PROCESO DE CAMINATA ALEATORIA

16100 C* L I T I L I Z A D A E N L A S I M U L A C I O N

16200

16300

16400 REAL INCR( 100)

16500

16600 C* I N I C I A L I Z A C I O N

\

106.

16700 J = O 16800 SANT = O

16900 s = o 17000 H = SQRT(DCAt1)

17 100 17200 C * 1N IC IP .N ITERACIONES

17300 DO 100 I I N C = 1 ,NINC 17400 J = J t l

17500 T I = IINC*DT

17600 DO 70 I = I N T ( ( T I - DT)/DCAM)+l ,INT(TI/DCAM)

17700 IF (PANCOM(SEM) .LE. O. 5 ) GOTO 60 17800 S = S + H

17900 GOTO 70

18000 60 S = S - H

18100 70 CONTINUE

18200 C" SE GUARDAN LOS INCREMENTOS

18300 INCR(J') = S - SANT

18400 SANT = S

18500 100 CONTINUE

18600 RETURN

18700 END 18800

18900

107.

PROGRAP.1AS

ESTO/GENERADOR

Este programa ut i l iza las subrut inas CAMINA y PTU, para generar y

g r a b a r en disco el arreglo de incrementos de MB de longitud DT:

( I N C ( J , K ) ; J = l , . . . , N I N C ; K = 1, ..., NTRAY)

donde N I N C = num. de incrementos NTRAY = n u m . de trayectorias

Las semi 11 as requeridas en 1 as trayectorias no se tienen que proporcionar, s i n o que se definen mediante la instrucción T I M E ( 1 ) .

19000

19 100

19200 C* PROGRAMA ESTO/GENERADOR

19300

19400 $RESET FREE 19500 .......................................................................

19600 C* ESTE PROGRANA OBTIENE N TRAYECTORIAS DE LONGITUD M DE MOVIMIENTO * 19700 C* BROWNIAN0 APROXIMADO POR CAMINATAS ALEATORIAS ESCALADAS, O POR * 19800 C* PROCESOS DE TRANSPORTE UNIFORPIE * i9900 C * S E U T I L I Z A N L A S S U B R U T I N A S : CAMINA(SEM,NINC,INCR,DT,DCAM) * 20000 c* PTU(SEM,NINC,INCR,DT,N) * 20100 C * Y GRABA EL ARREGLO DE INCREMENTOS * 20200 c* I N C ( J , K ) ; J=1, ..., M; K=1 , ..., N * 20300 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20400 F I L E 8 ( K I N D = D I S K , N E W F I L E )

20500

20600 R E A L I N C ( 100,500) ,MB( 100,500) , INCR( 100)

20700 DIMENSION ARCH( 6 )

103.

20800

20900 C* LECTURA DE PARAMETROS

2 1000 WRITE(6 ,/)"CUAL GENERADOR?? CA=O PTU=l"

2 1100 READ( 5 ,/) KGEN

2 1200 WRITE(G,/)"DAR NINC,TO,DT : I '

2 14QO WRITE(6,/)"DAR NUM. DE TRAYECTORIAS :"

2 1300 READ(5,/)NINCYTO,DT

2 1500 READ( 5 ,/)NTRAY 21600 35 CONTINUE 2 1700 IF(KGEN .EQ. o ) GOTO 40 2 1800 WRITE(G,/)"DAR VELOCIDAD DEL PTU N=? I'

2 1900 READ( 5 ,/)N 22000 GOTO 45

22100 40 WRITE(6,920) 2 2200 READ( 5 ,/) NSAL

22300 DCAM = DT/FLOAT( NSAL)

22400 45 CONTINUE

22500

22600 C* I N I C I A N ITERACIONES

22700 DO 100 K = l ,NTRAY

22800 SEM = TIME(1)

22900 IF(KGEN .EQ. 1) GOTO 60

23000 CALL CAMINA(SEM,NINC,INCR,DT,DCAPl) 23 100 GOTO 70

23200 60 CALL PTU(SEM,NINC,INCR,DT,N)

23300 70 CONTINUE

23400 C* SE FORMA LA FIATRIZ DE INCREMENTOS

2 3500 DO 80 J = l , N I N C

23600 I N C ( J ,K) = I N C R ( J )

23700 80 CONTINUE

23800 100 CONTINUE

23900 M = J 24000 C*SIGUE GRABACION O IMPRESION DE ARREGLOS

24100

109.

24200 C* GRABRCION DE LOS INCREMENTOS

24300 WRITE(6 ,/)"QUIERE GRABAR INCREMENTOS? NO=O, S I = 1 : ' I

24400 READ( 5 ,/) I N C G R A

24500 IF ( I N C G R A .EQ. O ) GOTO 110

24600 VEL = FLOAT(N)

2 4700 iF(KGEN .EQ. O ) VEL = NSAL

24800 WRITE( 6 ,/)"DAR A R C H I V O DE GRABACION DE INCREMENTOS : 'I

24900 REA~(5,910)(ARCH(Ij,I=ly5)

25000 ARCH(6) = I .

25 100 OPEN( 8,TITLE=ARCH ,KIND=DISK,FILETYPE=7)

25200 WRITE(8,/)NTRAY ,NINC,DT,VEL

25300 DO 105 J = l ,M 25400 WRITE( 8 , / ) ( INC( J ,K) ,K= l ,NTRAY

25500 105 CONTINUE

25600 LOCK( 8 )

25700 110 CONTINUE

25800 WRITE( 6,950)

25900 READ( 5 ,/) KOTRO

26000 IF(KOTR0 .EQ. 1 ) GOTO 35

26 100

26200 26300 C* SIGUEN FORMATOS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26400 910 FORMAT(6AG)

26500 920 FORMAT( ~ X , ~ D A R i: DE SALTOS POR INTERVALO DE MUESTREO : I )

26600 950 FORMAT(lX,'QUIERE OTRO PROCESO CON EL MISMO PARAM?? SI=l ,NO=O')

26700 END

26800

26900 C* SIGUEN LAS SUBRUTINAS PTU Y CAMINA ******************* 27000

27100

110.

ESTO/ESTADISTICAS/l

Este programa u t i l i za l as subru t inas CAMINA y PTU. Lee u n archivo en disco que contiene el arreglo

{INC(J,K); J=1, ... ,NINC; K = l , ... , N T R A Y )

generado y grabado anteriormente con el programa E S T O / G E N E R A D O R . Constru- ye NINC muestras de tamaño NTRAY de incrementos de MB, para cada una de las cuales obtiene e imprime de manera opcional:

l. Hi 2 tograma

2 . Media

3. Vari an za

4. Coef

5 . Coef

ic iente de asirnetria (sesgo)

ic ientes de Kurtosis

6 . Estadística de D'Agostino

7. Estadística de D'Agostino-Pearson.

Tznto la Kurtosis , como l a s dos últ imas estadísticas se uti l izan para probar la hipótesis de que l a s muestras de incrementos provienen de ur,a población normal , con media cero y vari anza DT.

Los parámetros que se requieren -y que el programa "pide" ai usua- r i o para poder obtener l as dos úl timas estadís t i cas , son DELTA y l/LAMBDA,

q u 2 aparecen en 1 a tab1 a de ID' Agostino , R. y E .S. Pearson ( 1 9 7 3 ) l , misma tab la que aparece como apéndice en [Seber, G.A.F. (op . c i t . ) ] . Para probar la hipótesis de normalidad mediante el coeficiente de Kurtosis, obtenemos

111.

la región de no rechazo de la f igura 1 de [D l Agostino , R. y E .S. Pearson (op. c i t . )I. La regi6n de no rechazo p a r a l a prueba de normalidad de D'Agostino, puede obtenerse a par t i r de los puntos percentiles del es ta- d í s t ico y cuya tabla aparece en [Seber, G .A.F. (op . c i t . )].

27200

27300

27400 C* PROGRAMA E S T O / E S T A D I S T I C A S / l

27500 27600 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27700 C* ESTE PROGRAMA LEE UN ARCHIVO DE N TRAYECTORIAS DE NMUEST PUNTOS * 27800 C* MANEJA EL ARCHIVO COMO NPWEST MUESTRAS DE LONGITUD N C/U * 27900 C* OBTIENE PARA CADA MUESTRA : * 28000 C* H ISTOGRAMA * 28100 C* MEDIA, VARIANZA, SESGO, KURTOSIS * 28200 C* Y ADEMAS E S T A D I S T I Z A S P A R A P R U E B A S DE NORNALIDAD DE: * 28300 C* D ' AGOSTINO * 28400 C* D 'AGOSTINO-PEARSON * 28500 C* OMNIBUS * 28600 C * E L A R C H I V O DE INCREMENTOS QUE SE LEERA, DEBE GRABARSE PREVIA.- * 28700 C' MENTE CON E L PROGRAMA ESTO/GENERADOR, Y ESTARA EN LA FORMA: * 28800 C * I I N C ( J , K ) ; J=1, ..., M; K=1, ..., N ) * 28900 ......................................................................

29000 F I L E

29100 F I L E

29200 29300

29400

29500

29600

29700 10 29800

29900

30000

G(KIND=REMOTE,MYUSE=IO,MAXRECSIZE=22) 7( KIND=PRINTER,FIAXRECSIZE=22)

D I M E N S I O N A ( 100,500) ,ARCH(6) ,FR( 100)

DIMENSION CAR( 120)

REAL RAN ,LONG ,MX ,MN ,KURT0

REAL RESTA ( 2 0 , 10 ) DATA CAR/120* ' * I / CONTINUE

W R I T E ( 6 , 9 1 1 )

R E A D ( 6 , 9 1 2 ) ( A R C H ( J ) ,J=1,5) A R C H ( 6 ) = ' . '

112.

30 100 WRITE(6,905)

30200 READ( 6 ,/ j KOMN I

30 300 IF(KOMN1 .EQ. O ) GOTO 70

30400 WRITE( 6,913)

30500 READ(6,/)ZDELYZLAM

30600 70 CONTINUE

30 700 WRITE( 6,915)

30800 READ(6 ,/) KGRA

30900 I F ( KGRA .EQ. O ) GOTO 77

31000 WRITE(6,701)

31100 READ(G,/)NI

31200 WRITE( 7,912) (ARCH(J) ,J=1,5)

31300 I F ( N I .EQ. 0 ) GOTO 6 0 1

31400 GOTO 77

31500 WRITE(G,/)"DAR EL NUMERO DE MUESTRAS"

31600 READ(6 ,/)NMUEST

31700 WRITE(G,/)"DAR EL NUMERO DE TRAYECTORIAS : I '

31800 READ(6 , / )N

31900 WRITE( 6 ,/)"DAR VELOCIDAD : I'

32000 READ(G,/)VEL

32100 77 CONTINUE

32200 OPEN(8,TITLE=ARCH ,KIND=DISK,FILETYPE=7)

32300 READ(8,/)N,NMUESTyDT,VEL 32 400 DO 501 J=l,NMUEST

32500 REP4D(8,/)(A(K,J) ,K=l,N)

32600 501 CONTINUE

32 700 LOCK( 8 )

32800 C* ELECCION DE LA UNIDAD DE SALIDA

32900 WRITE(6 ,/)"UNIDAD DE SALIDA?? TERMINAL=6 , IMPRESORA=7"

33000 READ( 6 ,/) I S

33100 C* IMPRESION DE ENCABEZADO

33200 WRITE( I S ,931)

33300 WRITE( I S ,930) (ARCH( I ) , I = l , 6 )

33400 WRITE( I S ,932)NMUEST

33500 WRITE( I S ,933)N

113.

33600 WRITE( I S ,935) DT

33700 WRITE( I S ,934)VEL

33800 C* ELECCION DE ALTURA DE? HISTOGRAMA

33900 ALTURA = 20

34000 I F ( I S .EQ. 6 ) GOTO 110

34100 IF(KGRA .EQ. O ) GOTO 110

34200 WRITE ( 6 ,/) "DAR ALTURA DEL H ISTOGRAMA ((6 1) : I '

34300 READ( 6 ,/)ALTURA

34400 110 CONTINUE

34500 C* I N I C I A N ITERACIONES POR MUESTRAS

34600

34 700

34800

34900

35000

35100

35200

35300

35400

35500

35600

35700

35800

35900

36000

36100

36200

36300

36400

36500

36600

DO 502 K=l,NMUEST

IF(KGRA .EQ. 1 ) WRITE(IS,~O~)K

IF(KGRA .EQ. o ) GOTO 350

DO 506 I=l,NI

FR( I)=O.O

506 CONTINUE

MX=-500 ;MN=500

DO 720 1=1 ,N I F ( M X .LT. A(1,K)) MX=A(I,K)

I F (MN .GT. A( 1 ,K)) MN=A( I ,K)

720 CONTINUE

RAN=MX-MN

LONG=RAN/NI

DO 725 J=l,N

I F (MN .LE. A(J,K).AND.MN+LoNG.GE.A(J,K))FR(l)=FR(l)+l 725 CONTINUE

DO 730 1=2 , N I

DO 740 J=l,N

I F (MN+( I- l )*LONG.LT.A(J ,K) .AND.MN+I*LONG.GE.A(J , K ) ) F R ( I ) = F R ( I ) + 1

740 CONTINUE

730 CONTINUE

36700 C* BUSQUEDA DE FREC. MAX. PARA ESCALAR HISTOGRAMA 36800 FMAX = O

36900 DO 220 I = l , N I

114

37000 IF(FMAX.LT.FR(I)) FMAX = FR(I)

37100 220 CONTINUE

37200 C* S I FMAX < ALTURA ENTONCES DEJAMOS EL HISTOGRAMA CON A L T U R A O R I G I N A L

37300 I F ( ALTURA. GE . FMAX) ALTURA = FMAX

37400 DO 750 I = 1 , N I

3 7500 GRA = F R ( I ) *ALTURA/FMAX

3 7600 JGRA = I N T ( G R A ) + 1 37700 J = F R ( I) 37800 I F ( J . E Q . 0 ) GOTO 800

37900 W R I T E ( I S ,760) I ,FR( I ) ,MN+( I - l ) * L O N G , ( C A R ( K K ) , K K = l , J G R A )

38000 GOTO 750

38100 800 W R I T E ( IS ,761.) I ,J ,MN+( I - l ) *LONG,MN+I *LONG

38200 750 CONTINUE

38300 C* S IGUE CALCULO DE ESTADISTICOS ************ 38400 350 CONTINUE

38500 I F ( KGRA. EQ. 1) W R I T E ( IS ,918) 38600

38700 C* CALCULO DE MEDIA, VARIANZA, SESGO Y KURTOSIS

38800 S1 = o; S2 = o; S3 = o; S4 = o 38900 DO 400 J=1 ,N

39000 S 1 = S1 + A ( J ,K)

39 100 S2 = S 2 + A ( J , K ) * * 2

39200 400 CONTINUE

39 300 PROM = S l / N

39400 VAR = (S2 - N * P R O M * * Z ) / ( N - l )

39500 IF( KGRA.EQ. 1) WRITE( IS , ~ ~ O ) P R O M , V A R

39600 C* CALCULO DE SESGO Y K U R T O S I S

39 700 DO 425 J = l , N

39800 S3 = S3 + ( A ( J , K ) - P R O M ) * * 3

39900 S4 = S 4 + (A (J ,K ) -PROM)* *4

40000 425 CONTINUE 40 100 DEN0 = VAR* (N-1 )

40200 SESGO = SQRT( N ) *S3/ ( DENO** l .5)

40300 KURT0 = F L O A T ( N ) *S4/ (DENO**2)

115.

40400 IF(KGRA.EQ.~) WRITE(IS,~~~)SESGO,KURTO 40500

40600 C* ESTADISTICOS PARA PRUEBAS DE NORMALIDAD ***************** 40700

40800 C* PRUEBA DE D' AGOSTINO. SEBER (1984) PAG. 144 40900 C* REORDENACION DE LA MUESTRA 41000 DO 440 LS=1, N-1 41 100 DO 435 J=2, N-LS+l 41200 IF(A(J-1 ,K) .LT.A(J ,K)) GOTO 435 41300 AUX = A(J-1,K) 41400 A(J-I,K) = A(J,K) 41500 A(J ,K) = AUX 41600 435 CONTINUE 41700 440 CONTINUE 41800 C* CALCULO DE D 41900 D = O

42000 DO 450 J=l,N 42 100 D = D + (J-(N+l)/Z)*A(J,K) 42200 450 CONTINUE 42300 D = D/SQRT( DENO*N**3) 42400 C* Y ES LA VERSION ESTANDARIZADA DE D 42500 Y = SQRT(N)*(D-l/(2*SQRT(3.14159)))/0.02998598

42600 IF(KGRA.EQ.1) WRITE( IS,925)Y 42 700

42800 C* SIGUEN PRUEBAS OMNIBUS 42900 IF(KOMNI.EQ.O) GOTO 490 43000 C* PRUEBA DE D' AGOSTINO-PEARSON. SEBER (1984) ,PAG. 144 43 100 XDP = DAP (SESGO ,ZDEL ,ZLAM) 43200 IF( KGRA. EQ. 1 ) WRITE ( IS ,940) XDP 43300 940 FORMAT(lX,'D AGOSTINO-PEARSON **Z = ',F12.6)

43400

43500 C* PRUEBA OMNIBUS DE BOWMAN-SHENTON. SEBER( 1984) ,PAG. 144 43600

43700 GOTO 490

116.

43800

43900 O M N I l = DAP(SESGO,ZDEL,ZLAM)

44000 OMNI2 = DAP( KURTO ,ZDEL ,ZLAtl)

44100 OMNI = OMNI 1**2 + OMN I2**2

44200 WRITE( I S ,942)OMNI

44300 942 FORMAT( l X , ' O M N I B U S ********* JI(2) = I ,F12.6

44400

44500 C* TERMINAN PRUEBAS OMNIBUS ******* I

44600

44700 490 CONTINUE

44800 C* ARREGLO DE ESTADISTICOS RESULTANTES (RESTA)

44900 RESTA( K,1) = PROM

45000 RESTA(K,2) = VAR 45 100 RESTA( K,3) = SESGO

45200 RESTA( K,4) = KURTO

45300 RESTA( K,5) = Y 45 400 RESTA(K,6) = XDP

45500

45600 C* TEP.MINAN ESTADISTICOS

45700

45800 502 CONTINUE

45900 WRITE( IS ,949)

46000 WRITE( I S ,950)

46 100 DO 550 K = l ,NFIUEST

46200 WRITE(IS,952) K,(RESTA(K,IE) , IE=1,6)

46300 550 CONTINUE

46400

46500 WRITE(6,/)"QUIERE DAR OTRO ARCHIVO?? NO=O S I = 1 :"

46600 READ(6 ,/)OTRO

46700 IF(OTRO.EQ.1) GOT0 10

46800 601 CALL EXIT

46900 503 FORMAT( // ,2OX , 'MUESTRA # ' ,14,//)

47000 701 FORMAT( lX, 'NUMERO DE INTERVALOS; (O=NO) ' . / ) 47100 761 FORMAT(1X,I3,3X,I4,2XY2(F8.5,2x))

47200 760 FORMAT(1X,I3,3X,I4,2XY2(F8.5,2X) ,120Al)

117.

47300

47400

47500

47600

47 700

47800

47900

48000

48100

48200

48,100

48400

48500

48600

48700

48800

48900

49000

911 FORMAT( l X , ' A R C H I V O : I ) 912 FORMAT( 6A6)

905 FORMAT( l X , ' Q U I E R E PRUEBA OMNIBUS?? NO=O, S I = 1 : )

913 FORNAT(lX,'DAR DELTA,l/LAMBDA (TABLA X EN SEBER) : I )

915 FORMAT(lX,'QUIERE HISTOGRAMAS?? NO=O, S I = 1 : I )

918 FORMAT(1XY'E5TADISTICOS : I )

920 FORMAT(lX,'MEDIA = ' ,F12.6,2XY'VARIANZA = ,F12.6)

922 FORMAT(lX,'SESCO = ' ,F12.6,2X,'KURTOSIS = ' ,F12.6)

925 FORMAT(lX,'PRUEBA DE NORMALIDAD DE DAGOSTINO ** Y = ' ,F12.6)

930 FORMAT( l X , ' A R C H I V O : ' :6A6)

931 FORMAT(//)

932 FORMAT( 1 X ,'NUM. DE MUESTRAS = ' , I 4 )

933 FORMAT( l X , ' N U M . DE TRAYECTORIAS = I , I 7 )

934 FORMAT( 1 X ,'VELOCIDAD = ' ,F8.1)

935 FORMAT(lX,'LONG. DE INTERVALO DE MlJESTREO = ,F12.6) 949 FORMAT( / / , 1QX , I ************ ESTADISTICAS k************** I

Y / )

950 FORMAT( 1X ,"MUEST. MEDIA VAR. SESGO KURTOSIS" , *I' Y ( D A) D-P" ,/)

49100 952

49200

49300 C* ** 49400

49500

49600

49700

49 800

49900

5000

FUNCTION DAP( Z ,ZDEL ,ZLAM)

ZY = Z*ZLAM

DAP = ZDEL*ALOG( ZY+SQRT( ZY*ZY+l . O ) )

RETURN

END

113.

ESTO/ESTADISTICAS/2

U t i l i z a n d o l a s s u b r u t i n a s CAMINA y PTU, es te p rog rama s imu la N t r a -

y e c t o r i as de MB , que p a r t e n de a+x, y c a l c u l a l a p r o p o r c i ó n de t r a y e c t o -

r i a s que se absorben por las barreras {X=a, X=b) ¶ p a r a

x = i ( b - a ) / N = i(RANGO)/N, i=l,.. . ,n-1

F ina lmente impr ime I X( I) , P( I ) , L I N F ( I ) , LSUP( I) , donde:

LINF(I) = X - 1 . 9 6 J m

LSUP(1) = x + 1.964-

Cuando l a p r u e b a s e r e a l i z a con un n i v e l de s i g n i f i c a n c i a i g u a l a

0.05 e l i n t e r v a l o ( L I N F ( I ) , LSUP( I ) ) es l a r e g i ó n de no rechazo de l a

h i p j t e s i s s i g u i e n t e :

La p ropo rc ión de t r a y e c t o r i a s que p a r t e n de a+x y que se absorben

p o r l a b a r r e r a (X=b) es i g u a l a x.

119.

5 1400

5 1500

51600 FILE 7 ( KIND=?RINTER,MAXRECSIZE=ZZ) 5 1700 DIMENSION TPSPOS( 200) ,TPSNEG(200) ,BAPOS(200) ,BANEG( 200)

5 180C DIMENS I O N ARCH ( 3 )

51900 INTEGER VEL,FINF(200) ,FSUP(200)

52000 REAL MP , INCR( 10 )

52100

52200 C INGRESA U N I D A D DE SALIDA

52300 WRITE(6,915)

52400 READ( 5 J) I S

52500

52600 C LECTURA DE PARAMETROS DE SIMULACION

52 700 ARCH( 1) = ' C A . ' 52800 ARCH(2) ='PTU. ' 52900 WRITE(6,920)

53000 READ( 5 ,/)NTRAY ,DT

53100 WRITE ( 6 , 922)

53200 READ( 5 ,/) KGEN

53300 IF ( KGEN . EQ. 1 ) GOTO 100

53400 WRITE(6,926)

53500 READ( 5 ,/) NSAL

53600 DCAM = DT/FLOAT( NSAL)

53700 GOTO 120

53800 100 WRITE( 6,924)

53900 READ( 5 ,/) VEL

54000 120 CONTINUE

54100 C LECTURA DE PARAMETROS SOBRE BARRERAS

54200

54300 WRITE( 6,930)

54400 READ( 5 y / ) NB 54500 WRITE(6,932) 54600 READ( 5 J) RANGO

54700 C ASIGNACION DE BARRERAS 54800

120.

54900 DO 130 1=1 ,NB

55000 BAPOS( I) = I*RANGO/NB

55 100 BANEG( I ) = -BAPOS ( I ) 55200 130 CONTINUE

55300

55400 C EMPIEZAN TRAYECTORIAS

55500

55600 DO 777 K=l,NTRAY

55700 C INICIALIZACION DE TRAYECTORIA

55 800 SEN = TIME(1)

55900 I P O S = 1 56000 INEG = 1

56 100 MB = O

56200 J = O

56300 DO 300 I=l,NB

56400 TPSPOS(1) = -1. 56500 TPSNEG(1) = -1. 56600 300 CONTINUE

56700

56800 C INICIAN ITERACIONES POR PUNTOS MUESTRALES

56900 200 CONTINUE

5 7000 J = J + l

57100 IF(KGEN.EQ.1) GOTO 210 5 7200 CALL CAMINA(SEM,~,INCR,DT,DCAM) 57300 GOTO 215 57400 210 CALL PTU(SEM,I,INCR,DT,VEL)

57500 215 MB = MB + INCR(1)

57600

57700 230 CONTINUE

5 7800 IF(MB.LT.BAPOS(IPOS)) GOTO 240 57900 TPSPOS( I P O S ) = J*DT

58000 IPOS = I P O S + 1 58100 IF ( IPOS.GT.NB) GOTO 250

58200 GOTO 230

58300 240 CONTINUE

121.

5 8400 IF(MB.GT.BAENG( INEG)) GOTO 200

58500 TPSNEG( INEG) = J*DT

58600 I N E G = I N E G + 1

58700 IF(INEG.LE.NB) GOTO 230

58800 250 CONTINUE

58900 C DETERMINACION DE BARRERAS DE SALIDA

59000

59100

59200 DO 320 I = l , N B - 1

59300 I B A P = I

59400 I B A N = NB - I

59500 IF(TPSPOS( IBAP) .LT.O) GOTO 305

59600 IF(TPSNEG( IBAN) .LT.O) GOTO 310

59700 IF(TPSPOS( IBAP) .LT.TPSNEG( IBAN)) GOTO 310

59800 305 FINF(1) = F INF(1 ) + i 59900 GOTO 320

60000 310 FSUP( I ) = FSUP( I ) + 1

60100 320 CONTINUE

60200 777 CONTINUE

60300

60400 C IMPRESION DE FRECUENCIAS DE SALIDAS POR BARRERAS

60500 C SUPERIOR E INFERIOR

60600 IF(KGEN.EQ.0) VEL = NSAL

60700 WRITE( I S ,950)

60800 WRITE(IS,~~~) ARCH(KGEN+I)

60900 WRITE( I S ,954) NTRAY ,DT ,VEL

61000 WRITE('IS,956)NB ,RANGO

6 1100 WRITE( I S ,960)

6 1200 DO 400 I=l ,NB- l

61300 X = BA?OS( I )

6 1400 DE = SQRT(X"(l.0-X)/FLOAT(NTRAY)) 6 1500 OINF = X-1.96*DE 6 1600 OSUP = X+1.96*DE

6 1700 P R I = FLOAT(FI!!F( I))/FLOAT(NTRAY)

122.

6 1800 PRS = F L O A T ( F S U P ( I ) ) / F L O A T ( N T R A Y )

6 1900 W R I T E ( I S ,965) I ,X ,PRI ,OINF,OSUP

62000 400 CONTINUE

62100 C* **** SIGUEN FORMATOS **** 62200

62300 915 F O R M A T ( l X , ' D A R U N I D A D DE S A L I D A T E R M I N A L = 6 , I M P R E S O R A = 7 ' )

62400 920 FORMAT( lX , 'NTRAY,DT ? ? : ' ) 6 2 5 0 0 9 2 2 FORMAT( lX , 'CUAL GENERADOR?? CA=O, PTU=l : I )

62600 926 FORMAT( 1X , 'NUM. DE SALTOS POR I N T E R V A L O DE MUESTREO ? ? : ' ) 62700 924 FORMAT( lX , 'VELOCIDAD DEL PTU.?? : ' ) 62800 930 FORMAT( 1X , 'DAR NUM. DE BARRERAS S I M E T R I C A S E N CERO ? ? : ' ) 62900 932 FORMAT( lX, ' DAR RANFO DE TRAYECTORIAS ? ? : ' ) 63000 950 FORMAT( lX , 'CONTEO DE SALIDAS POR BARRERAS SUP. E I N F . ' ,/, 63100 1 1 X , ' X O ( I ) : D I S T . DE L A BARRERA SUP. AL PUNTO DE P A R T I D A ' , / ,

63200 2 1 X , ' ( P R I ( I) : PROP. DE S A L I D A S POR L A BARRERA INFERIOR ' ,/, 63300 3 1 X , ' P R S ( I ) : PROP. DE S A L I D A S POR L A BARRERA SUPERIOR' ,/) 63400 952 FORMAT( lX , 'MB . S IMULADO POR : ' , A 6 )

63500 954 FORMAT( lX , 'NUM. DE TRAYECTORIAS : ' ,16,/,

63600 1 1 X ' INTERVALO DE MUESTREO : ' ,F8.5 ,/ , 63700 2 1 X , ' V E L O C I D A D : ' , I 5 )

63800 956 FORMAT( lX , 'NUM. DE B A R R E R A S S I M E T R I C A S E N CERO : ' ,14,/ , 63900 l l X , ' D I S T A N C I A E N T R E B A R R E R A S : ' ,F6.2)

64000 960 FORMAT(ZX,'I XO(I) PRI(I) LINF(I) L S U P ( I ) ~ )

64100 965 FORMAT(1X,I2,2X,F6.2,2X,F6.3,2(2X,F7.4)) 64200 END

64300

64400 C* S IGUEN SUBRUTINAS PTU Y C A M I N A

64500

64600

#

A P E N D I C E B

A q u í se presenta un programa para generar una trayectoria aproximada de una difusión.

Uti1 iza las rutinas CAMINA y PTU (apéndice A) pari; ger;erar trayector ias aproximadas de movimiento browniano. También ut i l iza la rut ina .RKF45 [Forsythe e t a l . (1977)J para integración numérica.

C * * * * * * * * * * * * * * * SIMULA UNA TRAYECTORIA DE UNA DIFUSION.

EXTERNAL F F

WQITE(4r/) "SEMILLA " READ(6r/) S E M

R E A 0 ( 6 < / ) N;? WRITE($//) SALIDA (6=TERMINAL, i'=IMPRESORA) "

READ(6,/) N S

'dRiTE(6,/) "OPCION (l=PTU, 2=CA) I'

E l = l . O E - 5 E2=l . O E - 5 N=l WRITE(NS/?10) TOP(N0P) TI-0. O T Z = T 1 Y ( 1 )=X0 IN=l NP=O

l o o IF(NP.GT.O.AND.NP.LT.IO) GO T O 2 0 0 IF(NOP.EQ.1) CALL PTU(SEMtIC,WW,Ht4) IF(NOP.EQ.2) CALL C A M I N A ( S E M I ~ O I W W ~ H / O . O ~ ) NP=O

200 NP=NP+l DW=WW (NP) /ti

124.

420 G O TO 500

4 3 0 G O TO 700

440 WRITE(6r920) IN/TIrY(1) IF(TI.EQ.T2) G O T O 500 G O T O 100

450 GO T O 700

4 6 0 G O T O 700

470 WRITE(6/920) INrTl r Y ( 1 1

G O T O 100 IF(Tl.EQ.TZ) G O TO 500

480 G O T O 700

5 0 0

7 0 0

T2=T1 +H IF(Y(l).LT.Xl+E?) G O T O 8 0 0 IF(Y(l).GT.XZ-EP) G O TO 8 0 0 G O TO 100

8 0 0 CALL E X I T

END

SUBROUTIqE FF(TrYrYP)

RETURN

E N D

125.

B I B L I O G R A F I A

BILLINGSLEY , P . (1968). "Convergence o f P r o b a b i 1 i ty Measures". John Wiley. New York.

B I R D , R.B. STEWART & E.N. LIGHTFOOT ( 1 9 6 0 ) . " T r a n s p o r t Phenomena. John Wiley & Sons, New York.

BPAUN, E . (1984) . "Movimiento Browniano" . Cienc ia - 35, 97-116.

CLPRK. J.M.C. ( 1 9 7 0 ) . " T h e r e p r e s e n t a t i o n o f f u n c t i o n a l s o f brownian mot ion as s tochas t i c i n teg ra l s " . Anna l s o f Ma th . S ta t . 41 , 1282 - 1295.

-

COURANT, R . , K.O. FRIEDRICHS & H. LOWY ( 1 9 2 8 ) . " b b e r d i e p a r t i e l l e d i f ferenzena1eichun;m der Phusik" . Math. Ann. 100, 32-74. Reimpreso (en i n g l é s ) e n IBM J. Res. Dev. - 11 ( 1 9 6 7 ) , 2 1 4 - 2 3 4 7

COURANT, R. 8. D. HILBERT (1966). "Methods o f M a t h e m a t i c a l P h i s i c s " , v o l . I. I n t e r s c i e n c e P u b l i s n e r s I n c . New York.

CROW, J. & M. KIMURA (1970) . "An i n t r o d u c t i o n t o P o p u l a t i o n Genetics Theory". Harper and Row, New York.

D'AGOSTINO, R . B . & E.S. PEARSON ( 1 9 7 3 ) . " T e s t s f o r d e p a r t u r e f r o m n o r m a l i t y . E m p i r i c a l r e s u l t s f o r t h e d i s t r i b u t i o n o f b, and f i l " B i o m e t r i k a 60, 613-622.

D I M A S I , G. & W . RUNGGALDIER (1980 ) . " A cont inuous t ime approx ima- t i o n t o o p t i m a l n o n l i n e a r f i l t e r i n g " . Pags. 153-160 de A r c h e t t i , F. & M. Cug iano , eds . (1980) , "Numer ica l techn iques fo r s tochas t i c sys tems" . Nor th H o l l a n d P u b l i s h i n g Co., Amsterdam. (1981) Appl. Math. Optim. - 7, 233-245.

DOOB, J .L . (1953) . "Stochast ic Processes" . John Wi ley & Sons. New York.

EINSTEIN, A . ( 1 9 5 6 ) . " I n v e s t i g a t i o n s on the Theory o f the Brownian Mct ion" . Dover . Inc. , New York.

FORSYTHE, G.E. 9. M.A. MOLER (1977). "Computer methods f o r mathemat ica l computa t ions" . Pren t i ce Ha l l , New York.

FREED", D. (1971) . "Brown ian Mot ion and D i f fus ion" . Ho lden Day, San Francisco.

FRIEDMAN, A. ( 1 9 7 5 ) . " S t o c h a s t i c D i f f e r e n t i a l E q u a t i o n s a n d A p p l i c a t i o n s " , Vol. I . Acadevic Press, New Yo rk '

GOROSTIZA, L. & R. GRIEGO (1980) . "Rate of convergence o f u n i f o r m t r a n s p o r t p r o c e s s e s t o b r o w n i a n m o t i o n a n d a p p l i c a t i o n t o s t o c h a s t i c i n t e g r a l s " . S t o c h a s t i c s , - 3, 291-303.

126.

GRIEGO, R . , n . C . HEATH & A . RUIZ-MONCAYO. ( 1 9 7 1 ) . "Almost sure convergence o f uniform t r a n s p o r t p r o c e s s e s t o brownian motion". A n n . Math. S t a t . - 42, 1129-1131.

HEDRICK, P .W . ( 1 9 8 4 ) . " G e n e t i c s o f p o p u l a t i o n s ' S c i e n c e Books I n t e r n a t i o n a l , B o s t o n .

HEPNANDEZ, D . B . (1981) . "Caminatas Aleator ias y Movimiento Browniano". Monografías del I n s t i t u t o de Matemáticas # 9 . U N A M , Cd. de Méx i co .

HERNANDEZ, D.B . ( 1 9 8 5 ) . " P r o c e s o s E s t o c á s t i c o s , con a p l i c a c i o n e s a l f i l t r a j e " . Rep. In t . Secc ión de Control del Departamento de Ingenier ía E l é c t r i c a del CINVESTAV, Cd. de México.

HERNANDEZ, D.B. & L . J . ALVAREZ (1985) . "E l método de Monte Carlo" . Informe Monográfico del Departamento de Matemáticas. UAM-Iztapalapa, Cd. de Méx i co .

JANSSEN, R . ( 1 9 8 4 ) . " D i s c r e t i z a t i o n o f t h e N i e n e r P r o c e s s i n Dif ference Methods f o r S t o c h a s t i c D i f f e r e n t i a l E q u a t i o n s " . S t o c h . P r o c . Appl. , 18, 361-369.

KAC, M . ( 1 9 5 1 ) . "On some connections between Probabil i ty Theory and D i f f e r e n t i a l and Integral Equat ions" , Proc . Second Berkeley Symposium on P r o b a b i l i t y and Mathematical S t a t i s t i c s , pp. 189-2 '15 . Univers i ty of C a l i f o r n i a P r e s s , B e r k e l e y .

K A C , M . ( 1 9 6 6 ) . "Can one hear the shape o f a drum?". American Mathematical Monthly, - 73 , 1 -23 .

KAKUTANI, S. ( 1 9 4 4 ) . "Two dimensional brownian motion and harmonic funct ions" , Proc . Imp. Acad. Tokyo, - 20 , 706-714 .

KIMURA, M . ( 1970) . "S tochas t i c Processes in Popula t ion Genet i cs " , en K.. Kojima (ed. ) "Mathematical topics in popula t ion genet i cs " . Spr inger Ver1 a g , B e r l i n .

KIMURA, M . & T . OTHA ( 1 9 7 1 ) . " T h e o r e t i c a l a s p e c t s o f p o p u l a t i o n g e n e t i c s " . P r i n c e t o n U n i v e r s i t y P r e s s , Princeton.

KNIGT, F . B . ( 1 9 8 1 ) . " E s s e n t i a l s o f Brownian Motion and Dif fusion" . AMS/Mathematical Surveys # 18, Providence , RI .

KUSHNER, H . J . (1977) . "Probabi 1 i t y methods for approximations in s t o c h a s t i c c o n t r o l and f o r e l l i p t i c e q u a t i o n s " . Academic P r e s s , New York.

KURTZ, T . ( 1 9 8 1 ) . "Approximation of population processes" . Regional Conference Series in Msthematics # 3 6 , SIAM, Phi lade lphia .

LEVI, E . ( 1 9 8 0 ) . " T e o r í a y Métodos de l a s Matemáticas Apl icadas'l . Ed. U N A M , Cd. de México.

LOEVE , M . ( 1977) . "Probabi l i t y Theory" , 4a . ed ic ión . Spr inger Verlag, New York.

MARUYAMA, G . (1955). "Continuous Markov processes and s t o c h a s t equations". Rend. Circ . Mat. Palermo - 4 , 4 8 - 9 0 .

McSHANE , E . J . ( 1 9 7 4 ) . " S t o c h a s t i c c a l c u l u s and s t o c h a s t i c mode Academic Press , New York.

NELSON, E . ( 1972) . "Dynanica l Theor ies o f Brownian Mction". P r i n c e t o n U n i v e r s i t y P r e s s , P r i n c e t o n , N . J .

127.

i c

1 S " *

R A O , M . ( 1 9 7 7 ) . "Brownian Motion and Class ical Potent ia ¡ Theory" . Lecture Notes Ser ies , No. 47 , Matematisk Inst i tut/Aarhus Univers i tet , Aarhus.

ROUGHARDEN, J. ( 1 9 8 0 ) . "Popul a t ion Genet i cs and Evolut ionary Eco1 ogy: An Introduct ion" . Mac Mi 1 1 an , Masachsset ts .

RUIZ MONCAYO, A . ( 1980) . " In t roducc ión a l a P r o b a b i l i d a d " . Fondo de Cultura Económica , Cd. de México.

RhELIN, W . ( 1 9 8 2 ) . " N u m e r i c a l t r e a t m e n t o f s t o c h a s t i c d i f f e r e n t i a l equations". SIAM J . Numer. Anal. - 19, 604-613.

SEBER, G.A.F. ( 1984) . "Mul t ivar ia te Observat ions" . John Wiley, New York.

STROOK, D . & S.R.S. VARADHAN (1979) . "Mul t id imens ional d i f fus ion processes" . Springer Verlag. New York.

V A R G A S , C . ( 1 9 8 5 ) . "Métodos numéricos para ecuaciones diferenciales" . Notas del IV Coloquio del Departamento de Matemáticas del CINVESTAV, Cd. de México .

VASQUEZ ABAD, F. (1984) . "S imulación de Procesos de D i f u s i ó n , con apl i cac iones en e l s e c t o r e n e r g é t i c o " . T e s i s de Maestr ía en E s t a d i s t i c a e Invest igación de Operaciones (IIMAS-UNAM) , Cd. de México.

WAX , N . ( 1954) . "Se lec ted papers on n o i s e and s t o c h a s t i c p r o c e s s e s " . Dover, Inc. New York.

WONG, E . & B . HAJEK ( 1 9 8 5 ) . " S t o c h a s t i c P r o c e s s e s i n E n g i n e e r i n g Systems". Springer Verlag, New York.

APROXIMACION DE PROCESOS DE DIFUSION148.206.53.84/tesiuami/reportesok/UAMR0823.pdf · Entonces, la...

Documents

Transcript of APROXIMACION DE PROCESOS DE DIFUSION148.206.53.84/tesiuami/reportesok/UAMR0823.pdf · Entonces, la...