#******************************************************************
# General Remarks:
#******************************************************************
#=================================================================================================
# This R source file contains code to estimate the invariant level crossing function. It also has
# a function that generates data from SDEs.
# The level crossing function uses matrix methods to do fast computations. If the time series
# and the number of specified levels is large, the resultant matrix will be huge,
# and may not fit in memory.

# The help file with example code is called 'helpinvariantSCTcodeweb.txt'. 

# This software does not come with any garantees.
# When using this software for scientific work, please cite:

# Wagenmakers, E.-J., Molenaar, P. C. M., Grasman, R. P. P. P., Hartelman, P. A. I.,
# & van der Maas, H. L. J. (2004). Tranformation Invariant Stochastic Catastrophe Theory.
# Manuscript submitted for publication.

#=================================================================================================
#========================================================================================

#******************************************************************
# R Code:
#******************************************************************

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
gendat = function(ndat = 1000, a = 0, b = 2, c = -1, dt = 0.1, plotit = T,
dvar1 = 0, dvar2 = 0, dvar2a = 2, dvar3 = 1, dvar3a = 1)
{
#=========================================================================================
# Generates simulated data from an SDE with specified drift and diffusion function
#=========================================================================================
	dat     = array();
	dat[1]  = rnorm(1, 0, sqrt(dt)); 

	for (j in 1:(ndat-1))
	{
		drift   = (a+b*dat[j]+c*dat[j]^3); #e.g., drift = -0.2*dat[j] gives OU process
		# default values give cusp process
		
		diffvar = dvar1*abs(dat[j]) + dvar2*(sin(dvar2a*dat[j])+1) + dvar3*dvar3a;
		#3 options for diffusion variance; edit the above line if you want to try others.
		
		error   = rnorm(1, 0, sqrt(diffvar*dt));
		dat[j+1]= dat[j] + drift*dt + error; # Cobb & Zacks (1985), Eq. 1.14
	}
	if (plotit)
		plot(dat, type="l", lwd = 1, main ="Simulated Data", ylab = "Value", xlab = "Time", cex.lab = 1.0, font.lab = 2, cex.axis = 0.9, bty = "n");
	
	invisible(dat); #same as "return", but without printing to console    
} 
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
levelcross = function(x, nlevels = length(x), dt = 0.1, plotit = F)
{
#==========================================================================
# Calculate number of upcrossings and downcrossings.
# cf. Hartelman (1997).
#==========================================================================
	n  = length(x);
	b  = seq(min(x), max(x), length = nlevels);
	y  = matrix(x, length(x), nlevels);
	y  = t(t(y) > b);
	y  = y[-1,] != y[-nrow(y),];
	lc = apply(y, 2, sum);
	# add Hartelman constant, cf. thesis p. 111, Eq. 4.5.13
	lc = lc / ( 2*(n-1)*sqrt(dt/(2*pi)) );

	if (plotit)
	{
		z = expression("#crossings"/(n-1)*sqrt(2*Delta/pi));		
		plot(b, lc, mgp = c(2.5, 1, 0), type="l", lwd = 2, main ="Level Crossing Function", ylab = z, xlab = "x", cex.lab = 1.0, font.lab = 2, cex.axis = 0.9, bty = "n");
		# mpg = c(3, 1, 0) is default. first = axis labels!; middle = tick labels  
	}	
	invisible(lc); #same as "return", but without printing to console	
}
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++