An example of the calculation and displaying of Ito-Gauss plots. It is a measure of the variability of linear trends
What is the variance of the collection of all the linear components,
- Variability of Millennial‐Scale Trends in the Geomagnetic Axial Dipole, Buffett et al. (2019)
- Decadal global temperature variability increases strongly with climate sensitivity, Nijsse et al. (2019)
Import libraries for calculations and getting data.
library(ItoGaussPlot)
library(rmatio)For this example, we will consider the Ornstein-Ulhenbeck process,
where
data <- read.mat('tX_g1D1e0.dat')
time <- data$tX[1,]
amplitude <- data$tX[2,]
dataVars <- list(gamma=1.0, D=1.0, dt=0.001)
dataVars$Nt <- ceiling(tail(time, n=1)/dataVars$dt)
plot(time, amplitude, type = 'l',
main='Ornstein-Uhlenbeck integration',
xlab='Time', ylab='Amplitude')
Calculating the variability of trends as a function of windowing length is done with the stdWindowSlopeCast() function. The time vector, amplitude vector, and window length(s) are arguments.
windows <- 10^seq(-1,2,len=20)
slopeStds <- stdWindowSlopeCast(time, amplitude, windows)For the Ornstein-Ulhenbeck process, a closed form expression of slope variability can be derived.
where
windowsTheory <- 10^seq(-1,2,len=50)
slopeStdsTheory <- OUanalytic(windowsTheory, dataVars$gamma, dataVars$D)
slopeStdsTheoryLimit <- OUanalytic(windowsTheory, dataVars$gamma, dataVars$D, "smallW")Below is a plot comparing the numerical calculations and the theoretical values.
plot(windowsTheory, slopeStdsTheory,
log='xy', type = 'l',
main='Standard deviation of slope by windowing',
xlab='Window length, w', ylab=expression('Std. of slope, σ'[b]))
lines(windowsTheory, slopeStdsTheoryLimit, lty=2)
points(windows, slopeStds, col="red", pch=16)
legend("topright", legend=c("Theory", "Small w limit", "Calculation"),
col=c("black", "black", "red"), lty=c(1,3,0), pch=c(0,0,16),
pt.cex=c(0,0,1), lwd=c(1.5,1.5,1), cex=0.8)
The package uses a different algorthm for unevenly spaces time-series.
data <- read.mat('tX_g1D1e0NU.dat')
timeNU <- data$tX[1,]
amplitudeNU <- data$tX[2,]
dataVarsNU <- list(gamma=1.0, D=1.0, dt=0.001)
dataVarsNU$Nt <- ceiling(tail(time, n=1)/dataVars$dt)
par(mfrow=c(1,2)) # set the plotting area into a 1*2 array
plot(timeNU, amplitudeNU, type = 'l',
main='Ornstein-Uhlenbeck integration',
xlab='Time', ylab='Amplitude')
hist(diff(timeNU),
main='Histogram of time-steps',
xlab='Time-step', ylab='Counts')
The stdWindowSlopeCast() function is used again, but theuneven timesteps in the time vector cause a different method to be called.
slopeStdsNU <- stdWindowSlopeCast(timeNU, amplitudeNU, windows)Below is a plot comparing the numerical calculations and the theoretical values.
plot(windowsTheory, slopeStdsTheory,
log='xy', type = 'l',
main='Standard deviation of slope by windowing',
xlab='Window length, w', ylab=expression('Std. of slope, σ'[b]))
lines(windowsTheory, slopeStdsTheoryLimit, lty=2)
points(windows, slopeStdsNU, col="red", pch=16)
legend("topright", legend=c("Theory", "Small w limit", "Calculation"),
col=c("black", "black", "red"), lty=c(1,3,0), pch=c(0,0,16),
pt.cex=c(0,0,1), lwd=c(1.5,1.5,1), cex=0.8)
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent