11options(prompt = " R> " ,scipen = 16 ,digits = 4 ,warning = FALSE , message = FALSE )
2- library(Sim.DiffProc )
2+ library(Sim.DiffProc )
33
44# # -------------------------------------------------------------------
55
@@ -41,6 +41,46 @@ print(dfptsde1d(fpt1d))
4141plot(dfptsde1d(fpt1d ),hist = TRUE ,nbins = " FD" ) # # histogramm
4242plot(dfptsde1d(fpt1d )) # # kernel density
4343
44+ # # -------------------------------------------------------------------
45+
46+ f <- expression( (1 - 0.5 * x ) )
47+ g <- expression( 1 )
48+ mod1d <- snssde1d(drift = f ,diffusion = g ,x0 = 2.1 ,M = 20 ,type = " str" )
49+
50+ St <- expression(2 * (1 - sinh(0.5 * t )) )
51+ fpt1d <- fptsde1d(mod1d , boundary = St )
52+ print(fpt1d )
53+ head(fpt1d $ fpt , n = 4 )
54+
55+
56+ mean(fpt1d )
57+ moment(fpt1d , center = TRUE , order = 2 ) # # variance
58+ Median(fpt1d )
59+ Mode(fpt1d )
60+ quantile(fpt1d )
61+ kurtosis(fpt1d )
62+ skewness(fpt1d )
63+ cv(fpt1d )
64+ min(fpt1d )
65+ max(fpt1d )
66+ moment(fpt1d , center = TRUE , order = 4 )
67+ moment(fpt1d , center = FALSE , order = 4 )
68+
69+ summary(fpt1d )
70+
71+
72+ plot(time(mod1d ),mod1d $ X [,1 ],type = " l" ,lty = 3 ,ylab = " X(t)" ,xlab = " time" ,axes = F )
73+ curve(2 * (1 - sinh(0.5 * x )),add = TRUE ,col = 2 )
74+ points(fpt1d $ fpt [1 ],2 * (1 - sinh(0.5 * fpt1d $ fpt [1 ])),pch = 19 ,col = 4 ,cex = 0.5 )
75+ lines(c(fpt1d $ fpt [1 ],fpt1d $ fpt [1 ]),c(0 ,2 * (1 - sinh(0.5 * fpt1d $ fpt [1 ]))),lty = 2 ,col = 4 )
76+ axis(1 , fpt1d $ fpt [1 ], bquote(tau [S(t )]== .(fpt1d $ fpt [1 ])),col = 4 ,col.ticks = 4 )
77+ legend(' topleft' ,col = c(1 ,2 ,4 ),lty = c(1 ,1 ,NA ),pch = c(NA ,NA ,19 ),legend = c(expression(X [t ]),expression(S(t )),expression(tau [S(t )])),cex = 0.8 ,bty = ' n' )
78+ box()
79+
80+ print(dfptsde1d(fpt1d ))
81+ plot(dfptsde1d(fpt1d ),hist = TRUE ,nbins = " FD" ) # # histogramm
82+ plot(dfptsde1d(fpt1d ))
83+
4484
4585# # -------------------------------------------------------------------
4686
@@ -93,6 +133,54 @@ plot(denJ,display="persp",main="Bivariate Density of F.P.T",xlab=expression(tau[
93133
94134# # -------------------------------------------------------------------
95135
136+ fx <- expression(5 * (- 1 - y )* x , 5 * (- 1 - x )* y )
137+ gx <- expression(0.5 * y ,0.5 * x )
138+ mod2d <- snssde2d(drift = fx ,diffusion = gx ,x0 = c(x = - 1 ,y = 1 ),M = 20 ,type = " str" )
139+
140+ St <- expression(sin(2 * pi * t ))
141+ fpt2d <- fptsde2d(mod2d , boundary = St )
142+ print(fpt2d )
143+ head(fpt2d $ fpt , n = 4 )
144+
145+ mean(fpt2d )
146+ moment(fpt2d , center = TRUE , order = 2 ) # # variance
147+ Median(fpt2d )
148+ Mode(fpt2d )
149+ quantile(fpt2d )
150+ kurtosis(fpt2d )
151+ skewness(fpt2d )
152+ cv(fpt2d )
153+ min(fpt2d )
154+ max(fpt2d )
155+ moment(fpt2d , center = TRUE , order = 4 )
156+ moment(fpt2d , center = FALSE , order = 4 )
157+
158+ summary(fpt2d )
159+
160+ plot(ts.union(mod2d $ X [,1 ],mod2d $ Y [,1 ]),col = 1 : 2 ,lty = 3 ,plot.type = " single" ,type = " l" ,ylab = " " ,xlab = " time" ,axes = F )
161+ curve(sin(2 * pi * x ),add = TRUE ,col = 3 )
162+ points(fpt2d $ fpt $ x [1 ],sin(2 * pi * fpt2d $ fpt $ x [1 ]),pch = 19 ,col = 4 ,cex = 0.5 )
163+ lines(c(fpt2d $ fpt $ x [1 ],fpt2d $ fpt $ x [1 ]),c(sin(2 * pi * fpt2d $ fpt $ x [1 ]),- 10 ),lty = 2 ,col = 4 )
164+ axis(1 , fpt2d $ fpt $ x [1 ], bquote(tau [X [S(t )]]== .(fpt2d $ fpt $ x [1 ])),col = 4 ,col.ticks = 4 )
165+ points(fpt2d $ fpt $ y [1 ],sin(2 * pi * fpt2d $ fpt $ y [1 ]),pch = 19 ,col = 5 ,cex = 0.5 )
166+ lines(c(fpt2d $ fpt $ y [1 ],fpt2d $ fpt $ y [1 ]),c(sin(2 * pi * fpt2d $ fpt $ y [1 ]),- 10 ),lty = 2 ,col = 5 )
167+ axis(1 , fpt2d $ fpt $ y [1 ], bquote(tau [Y [S(t )]]== .(fpt2d $ fpt $ y [1 ])),col = 5 ,col.ticks = 5 )
168+ legend(' topright' ,col = 1 : 5 ,lty = c(1 ,1 ,1 ,NA ,NA ),pch = c(NA ,NA ,NA ,19 ,19 ),legend = c(expression(X [t ]),expression(Y [t ]),expression(S(t )),expression(tau [X [S(t )]]),expression(tau [Y [S(t )]])),cex = 0.8 ,inset = .01 )
169+ box()
170+
171+ denM <- dfptsde2d(fpt2d , pdf = ' M' )
172+ print(denM )
173+ plot(denM )
174+ plot(denM , hist = TRUE )
175+
176+ denJ <- dfptsde2d(fpt2d , pdf = ' J' ,n = 100 )
177+ print(denJ )
178+ plot(denJ ,display = " contour" ,main = " Bivariate Density of F.P.T" ,xlab = expression(tau [x ]),ylab = expression(tau [y ]))
179+ plot(denJ ,display = " image" ,main = " Bivariate Density of F.P.T" ,xlab = expression(tau [x ]),ylab = expression(tau [y ]))
180+ plot(denJ ,display = " persp" ,main = " Bivariate Density of F.P.T" ,xlab = expression(tau [x ]),ylab = expression(tau [y ]))
181+
182+ # # -------------------------------------------------------------------
183+
96184fx <- expression(4 * (- 1 - x )* y , 4 * (1 - y )* x , 4 * (1 - z )* y )
97185gx <- rep(expression(0.2 ),3 )
98186mod3d <- snssde3d(drift = fx ,diffusion = gx ,x0 = c(x = 2 ,y = - 2 ,z = 0 ),M = 20 )
@@ -136,6 +224,55 @@ print(denM)
136224plot(denM )
137225plot(denM ,hist = TRUE )
138226
227+ denJ <- dfptsde3d(fpt3d ,pdf = " J" )
228+ print(denJ )
229+ plot(denJ ,display = " rgl" )
230+
231+ # # -------------------------------------------------------------------
232+
233+ fx <- expression(4 * (1 - x )* y , 4 * (1 - y )* x , 4 * (1 - z )* y )
234+ gx <- rep(expression(0.2 ),3 )
235+ mod3d <- snssde3d(drift = fx ,diffusion = gx ,x0 = c(x = - 2 ,y = 2 ,z = - 2 ),M = 20 ,type = " str" )
236+
237+ St <- expression(- 1.5 + 3 * t )
238+ fpt3d <- fptsde3d(mod3d , boundary = St )
239+ print(fpt3d )
240+ head(fpt3d $ fpt , n = 4 )
241+
242+ mean(fpt3d )
243+ moment(fpt3d , center = TRUE , order = 2 ) # # variance
244+ Median(fpt3d )
245+ Mode(fpt3d )
246+ quantile(fpt3d )
247+ kurtosis(fpt3d )
248+ skewness(fpt3d )
249+ cv(fpt3d )
250+ min(fpt3d )
251+ max(fpt3d )
252+ moment(fpt3d , center = TRUE , order = 4 )
253+ moment(fpt3d , center = FALSE , order = 4 )
254+
255+ summary(fpt3d )
256+
257+ plot(ts.union(mod3d $ X [,1 ],mod3d $ Y [,1 ],mod3d $ Z [,1 ]),col = 1 : 3 ,lty = 3 ,plot.type = " single" ,type = " l" ,ylab = " " ,xlab = " time" ,axes = F )
258+ curve(- 1.5 + 3 * x ,add = TRUE ,col = 4 )
259+ points(fpt3d $ fpt $ x [1 ],- 1.5 + 3 * fpt3d $ fpt $ x [1 ],pch = 19 ,col = 5 ,cex = 0.5 )
260+ lines(c(fpt3d $ fpt $ x [1 ],fpt3d $ fpt $ x [1 ]),c(- 1.5 + 3 * fpt3d $ fpt $ x [1 ],- 10 ),lty = 2 ,col = 5 )
261+ axis(1 , fpt3d $ fpt $ x [1 ], bquote(tau [X [S(t )]]== .(fpt3d $ fpt $ x [1 ])),col = 5 ,col.ticks = 5 )
262+ points(fpt3d $ fpt $ y [1 ],- 1.5 + 3 * fpt3d $ fpt $ y [1 ],pch = 19 ,col = 6 ,cex = 0.5 )
263+ lines(c(fpt3d $ fpt $ y [1 ],fpt3d $ fpt $ y [1 ]),c(- 1.5 + 3 * fpt3d $ fpt $ y [1 ],- 10 ),lty = 2 ,col = 6 )
264+ axis(1 , fpt3d $ fpt $ y [1 ], bquote(tau [Y [S(t )]]== .(fpt3d $ fpt $ y [1 ])),col = 6 ,col.ticks = 6 )
265+ points(fpt3d $ fpt $ z [1 ],- 1.5 + 3 * fpt3d $ fpt $ z [1 ],pch = 19 ,col = 7 ,cex = 0.5 )
266+ lines(c(fpt3d $ fpt $ z [1 ],fpt3d $ fpt $ z [1 ]),c(- 1.5 + 3 * fpt3d $ fpt $ z [1 ],- 10 ),lty = 2 ,col = 7 )
267+ axis(1 , fpt3d $ fpt $ z [1 ], bquote(tau [Z [S(t )]]== .(fpt3d $ fpt $ z [1 ])),col = 7 ,col.ticks = 7 )
268+ legend(' topright' ,col = 1 : 7 ,lty = c(1 ,1 ,1 ,1 ,NA ,NA ,NA ),pch = c(NA ,NA ,NA ,NA ,19 ,19 ,19 ),legend = c(expression(X [t ]),expression(Y [t ]),expression(Z [t ]),expression(S(t )),expression(tau [X [S(t )]]),expression(tau [Y [S(t )]]),expression(tau [Z [S(t )]])),cex = 0.8 ,inset = .01 )
269+ box()
270+
271+ denM <- dfptsde3d(fpt3d , pdf = " M" )
272+ print(denM )
273+ plot(denM )
274+ plot(denM ,hist = TRUE )
275+
139276denJ <- dfptsde3d(fpt3d ,pdf = " J" )
140277print(denJ )
141278plot(denJ ,display = " rgl" )
0 commit comments