R/CauchyDGP.R

In kavehsn/PredictiveDGP: Predictive DGPs exhibiting various distributions

#' Cauchy errors 
#' 
#' 
#' The processes \eqn{y} and \eqn{x} are generated by a predictive regression of the form \eqn{y_t=\beta x_{t-1}+\varepsilon_t} for \eqn{t=1,\cdots,T},
#' in which the regressors follow an AR(1) process - i.e. \eqn{x_t=\theta x_{t-1}+u_t}.  
#' The predictor's errors are distributed according to \eqn{u_t\sim N(0,1)}, whereas the disturbances 
#' of the predictive regression, \eqn{\varepsilon_t}, are distributed \eqn{\varepsilon_t\sim Cauchy}. An example of a predictive regression
#' DGP with Cauchy perturbations can be found in \insertCite{dufour1995exact;textual}{PredictiveDGP}.
#' The initial value of the process \eqn{x} is generated by \eqn{x_0=\frac{w_0}{\sqrt{1-\theta^2}}}, 
#' where \eqn{w_t\sim N(0,1)}. Finally, the contemporaneous correlation between the disturbances \eqn{\varepsilon_t} and \eqn{u_t} is captured by \eqn{\rho\varepsilon_t+w_t\sqrt{1-\rho^2}}.
#' 
#' 
#' @param n the number of observations
#' @param beta the regressor coefficient of the predictive regression
#' @param theta the autocorrelation coefficient of the predictor
#' @param rho the contemporaneous correlation coefficient    
#' @keywords endogeneity persistency Cauchy
#' @import pracma
#' @export 
#' @examples
#' CauchyDGP(n=50, beta=0.5, theta=0.999, rho=0.9)
#' @references
#' \insertAllCited{}
#' @importFrom Rdpack reprompt
#
CauchyDGP <- function(n, beta, theta, rho){

	#Generate the Cauchy error terms for the predictive regression, and the standard normal errors for the regressor
	eps <- rt(n, df=1)
	w <- rnorm(n, 0, 1)

	#Generate empty variables
	u <- rep(0, times=n)
	x <- rep(0, times=n)
	y <- rep(0, times=n)
	
	#Generate the initial predictor
	x[1] <- w[1]/(sqrt(1-theta^2))

	#Generate the error terms for the regressors. Note the possibility of feedback from the predictive regression errors to future x.
	for(i in 1:n){

		u[i] <- (rho*eps[i])+(w[i]*sqrt(1-rho^2))
	}
	
	#Generate the X and Y vectors
	for(k in 2:n){

		y[k] <- beta*x[k-1]+eps[k]
		x[k] <- theta*x[k-1]+u[k]

	}

	#Outcome variable
	z <- flipud(cbind(y, x))
	
	#Return output
	return(z)

}