R/CFAR.r

Defines functions g_cfar1 g_cfar2 g_cfar g_cfar2h est_cfar F_test_cfar F_test_cfarh est_cfarh p_cfar p_cfar_part

Documented in est_cfar est_cfarh F_test_cfar F_test_cfarh g_cfar g_cfar1 g_cfar2 g_cfar2h p_cfar p_cfar_part

#' Generate a CFAR(1) Process
#'
#' Generate a convolutional functional autoregressive process with order 1.
#' @param tmax length of time.
#' @param rho parameter for O-U process (noise process). 
#' @param phi_func convolutional function. Default is density function of normal distribution with mean 0 and standard deviation 0.1.
#' @param grid the number of grid points used to constrcut the functional time series. Default is 1000.
#' @param sigma the standard deviation of O-U process. Default is 1.
#' @references
#' Liu, X., Xiao, H., and Chen, R. (2016) Convolutional autoregressive models for functional time series. \emph{Journal of Econometrics}, 194, 263-282.
#' @return The function returns a list with components:
#' \item{cfar1}{a tmax-by-(grid+1) matrix following a CFAR(1) process.}
#' \item{epsilon}{the innovation at time tmax.}
#' @examples
#' phi_func= function(x)   
#' {
#'  	return(dnorm(x,mean=0,sd=0.1))
#' }
#' y=g_cfar1(1000,5,phi_func)
#' @export
g_cfar1 <- function(tmax=1001,rho=5,phi_func=NULL,grid=1000,sigma=1){
	#################################
	###Simulate CFAR(1) processes
	#############################
	###parameter setting
	#rho	### parameter for error process epsilon_t, O-U process
	#sigma is the standard deviation of OU-process
	#tmax	### length of time
	#iter	### number of replications in the simulation, iter=1
	#grid	### the number of grid points used to construct the functional time series X_t and epsilon_t
        # phi_func:   User supplied convolution function phi(.). The following default is used if not specified.
	if(is.null(phi_func)){
		phi_func= function(x){
			return(dnorm(x,mean=0,sd=0.1))
            }
    	}
	if(is.null(grid)){
		grid=1000
	}
	if(is.null(sigma)){
		sigma=1
	}
	#########################################################################################################
	#### OUTPUTS
	#######################
	#### A matrix f_grid is generated, with (tmax) rows and (grid+1) columns. 
	########################################################################################################

	#################################################################################################
	x_grid<- seq(0, 1, by= 1/grid); ### the grid points for input variables of functional time series in [0,1]
	b_grid<- seq(-1, 1, by=1/grid);	### the grid points for input variables of b-spline functions in [-1,1]
	eps_grid <- matrix(0, 1, grid+1);	### the grid points for input variables of error process in [0,1]
	fi <- exp(-rho/grid);		### the correlation of two adjacent points for an O-U process
	phi_b <- phi_func(b_grid)		### phi(s) when s in [-1,1]
	### the convolutional function values phi(s-u), for different s in [0,1]
	phi_matrix <- matrix(0, (grid+1), (grid+1));
	for (i in 1:(grid+1)){
		phi_matrix[i,]= phi_b[seq((i+grid),i,by=-1)]/(grid+1)
	}

	##############################
	### Generate data
	### functional time series f_grid is X_t in the paper
	#### A matrix f_grid is generated, with (tmax*iter) rows and (grid+1) columns. 
	#### It contains iter CFAR(1) processes. The i-th process is in rows (i-1)*tmax+1:tmax.

	f_grid=matrix(0,tmax,grid+1)
	i=1
	eps_grid= arima.sim(n= grid+1, model=list(ar=fi),rand.gen= function(n)sqrt((1-fi^2))*rnorm(n,0,1))  ###error process
	f_grid[1,]= eps_grid;	
	for (i in 2:tmax){
		eps_grid= arima.sim(n= grid+1, model=list(ar=fi),rand.gen= function(n)sqrt(1-fi^2)*rnorm(n,0,1))
		f_grid[i,]= phi_matrix%*%f_grid[i-1,]+ eps_grid;
	}
	g_cfar1 <- list(cfar1=f_grid*sigma,epsilon=eps_grid*sigma)
	return(g_cfar1)
}


#' Generate a CFAR(2) Process
#'
#' Generate a convolutional functional autoregressive process with order 2.
#' @param tmax length of time.
#' @param rho parameter for O-U process (noise process). 
#' @param phi_func1 the first convolutional function. Default is 0.5*x^2+0.5*x+0.13.
#' @param phi_func2 the second convolutional function. Default is 0.7*x^4-0.1*x^3-0.15*x.
#' @param grid the number of grid points used to constrcut the functional time series. Default is 1000.
#' @param sigma the standard deviation of O-U process. Default is 1.
#' @references
#' Liu, X., Xiao, H., and Chen, R. (2016) Convolutional autoregressive models for functional time series. \emph{Journal of Econometrics}, 194, 263-282.
#' @return The function returns a list with components:
#' \item{cfar2}{a tmax-by-(grid+1) matrix following a CFAR(1) process.}
#' \item{epsilon}{the innovation at time tmax.}
#' @examples
#' phi_func1= function(x){
#'	return(0.5*x^2+0.5*x+0.13)
#' }
#' phi_func2= function(x){
#'	return(0.7*x^4-0.1*x^3-0.15*x)
#' }
#' y=g_cfar2(1000,5,phi_func1,phi_func2)
#' @export
g_cfar2 <- function(tmax=1001,rho=5,phi_func1=NULL, phi_func2=NULL,grid=1000,sigma=1){
	#################################
	###Simulate CFAR(2) processes

	#########################################################################################
	#### OUTPUTS
	#######################
	#### A matrix f_grid is generated, with (tmax) rows and (grid+1) columns.
	###########################################################################################

	#############################
	###parameter setting
	## rho=5	### parameter for error process epsilon_t, O-U process
	##tmax= 1001;	### length of time
	##iter= 100;	### number of replications in the simulation
	##grid=1000;	### the number of grid points used to construct the functional time series X_t and epsilon_t
        ## phi_func1 and phi_func2: User specified convolution functions. The following defaults are used if not given
	if(is.null(phi_func1)){
		phi_func1= function(x){
			return(0.5*x^2+0.5*x+0.13)
    		}
	}
	if(is.null(phi_func2)){
	    phi_func2= function(x){
		  return(0.7*x^4-0.1*x^3-0.15*x)
	    }
	}
	if(is.null(grid)){
		grid=1000
	}
	if(is.null(sigma)){
		sigma=1
	}
	#######################################################################################################
	x_grid= seq(0, 1, by= 1/grid);	### the grid points for input variables of functional time series in [0,1]
	b_grid= seq(-1, 1, by=1/grid);	### the grid points for input variables of b-spline functions in [-1,1]
	eps_grid= matrix(0, 1, grid+1);	### the grid points for input variables of error process in [0,1]
	fi = exp(-rho/grid);		### the correlation of two adjacent points for an O-U process
	phi_b1= phi_func1(b_grid)	### phi_1(s) when s in [-1,1]
	phi_b2= phi_func2(b_grid)	### phi_2(s) when s in [-1,1]

	### the convolutional function values phi_1(s-u) and phi_2(s-u), for different s in [0,1]
	phi_matrix1= matrix(0, (grid+1), (grid+1));
	for (i in 1:(grid+1)){
		phi_matrix1[i,]= phi_b1[seq((i+grid),i,by=-1)]/(grid+1)
	}
	phi_matrix2= matrix(0, (grid+1), (grid+1));
	for (i in 1:(grid+1)){
		phi_matrix2[i,]= phi_b2[seq((i+grid),i,by=-1)]/(grid+1)
	}
	##############################
	### Generate data
	### functional time series f_grid is X_t in the paper
	#### A matrix f_grid is generated, with (tmax*iter) rows and (grid+1) columns. 
	##  It contains iter CFAR(2) processes. The i-th process is in rows (i-1)*tmax+1:tmax.

	f_grid=matrix(0,tmax,grid+1)
	i=1
	eps_grid= arima.sim(n= grid+1, model=list(ar=fi),rand.gen= function(n)sqrt(1-fi^2)*rnorm(n,0,1))
	f_grid[1,]= eps_grid;
	i=2
	eps_grid= arima.sim(n= grid+1, model=list(ar=fi),rand.gen= function(n)sqrt(1-fi^2)*rnorm(n,0,1))
	f_grid[i,]= phi_matrix1%*%as.matrix(f_grid[i-1,])+ eps_grid;
	for (i in 3:tmax){
		eps_grid= arima.sim(n= grid+1, model=list(ar=fi),rand.gen= function(n)sqrt(1-fi^2)*rnorm(n,0,1))
		f_grid[i,]= phi_matrix1%*%as.matrix(f_grid[i-1,])+ phi_matrix2%*%as.matrix(f_grid[i-2,])+ eps_grid;
	}
	### The last iteration of epsilon_t is returned.
	g_cfar2 <- list(cfar2=f_grid*sigma,epsilon=eps_grid*sigma)
}


#' Generate a CFAR Process
#'
#' Generate a convolutional functional autoregressive process.
#' @param tmax length of time.
#' @param rho parameter for O-U process (noise process). 
#' @param phi_list the convolutional function(s). Default is the density function of normal distribution with mean 0 and standard deviation 0.1.
#' @param grid the number of grid points used to constrcut the functional time series. Default is 1000.
#' @param sigma the standard deviation of O-U process. Default is 1.
#' @references
#' Liu, X., Xiao, H., and Chen, R. (2016) Convolutional autoregressive models for functional time series. \emph{Journal of Econometrics}, 194, 263-282.
#' @return The function returns a list with components:
#' \item{cfar}{a tmax-by-(grid+1) matrix following a CFAR(1) process.}
#' \item{epsilon}{the innovation at time tmax.}
#' @examples
#' phi_func= function(x)   
#' {
#'	return(dnorm(x,mean=0,sd=0.1))
#' }
#' test=g_cfar(1000,5,phi_func)
#'
#' phi_func1= function(x){
#'	return(0.5*x^2+0.5*x+0.13)
#' }
#' phi_func2= function(x){
#'	return(0.7*x^4-0.1*x^3-0.15*x)
#' }
#' phi_list=list(phi_func1,phi_func2)
#' y=g_cfar(1000,5,phi_list)
#' @export
g_cfar <- function(tmax=1001,rho=5,phi_list=NULL, grid=1000,sigma=1){
	#################################
	###Simulate CFAR processes

	#########################################################################################
	#### OUTPUTS
	#######################
	#### A matrix f_grid is generated, with (tmax) rows and (grid+1) columns.
	###########################################################################################

	#############################
	###parameter setting
	## rho=5	### parameter for error process epsilon_t, O-U process
	##tmax= 1001;	### length of time
	##iter= 100;	### number of replications in the simulation
	##grid=1000;	### the number of grid points used to construct the functional time series X_t and epsilon_t
       ## phi_func User specified convolution functions. The following defaults are used if not given
	if(is.null(phi_list)){
		phi_func1= function(x){
			return(dnorm(x,mean=0,sd=0.1))
    		}
		phi_list=phi_func1
	}
	if(is.null(grid)){
		grid=1000
	}
	if(is.null(sigma)){
		sigma=1
	}
	#######################################################################################################
	x_grid= seq(0, 1, by= 1/grid);	### the grid points for input variables of functional time series in [0,1]
	b_grid= seq(-1, 1, by=1/grid);	### the grid points for input variables of b-spline functions in [-1,1]
	eps_grid= matrix(0, 1, grid+1);	### the grid points for input variables of error process in [0,1]
	fi = exp(-rho/grid);		### the correlation of two adjacent points for an O-U process
	p=length(phi_list)
	if(is.list(phi_list)){
		phi_b=sapply(phi_list,mapply,b_grid)	
	}else{
		phi_b=matrix(phi_list(b_grid),2*grid+1,1)
	}
	### the convolutional function values phi_1(s-u) and phi_2(s-u), for different s in [0,1]
	phi_matrix= matrix(0, (grid+1), p*(grid+1));
	for(j in 1:p){
		for (i in 1:(grid+1)){
			phi_matrix[i,((j-1)*(grid+1)+1):(j*(grid+1))]= phi_b[seq((i+grid),i,by=-1),j]/(grid+1)
		}
	}

	##############################
	### Generate data
	### functional time series f_grid is X_t in the paper
	#### A matrix f_grid is generated, with (tmax*iter) rows and (grid+1) columns. 
	##  It contains iter CFAR(2) processes. The i-th process is in rows (i-1)*tmax+1:tmax.

	f_grid=matrix(0,tmax,grid+1)
	i=1
	eps_grid= arima.sim(n= grid+1, model=list(ar=fi),rand.gen= function(n)sqrt(1-fi^2)*rnorm(n,0,1))
	f_grid[1,]= eps_grid;
	if(p>1){
		for(i in 2:p){
			eps_grid= arima.sim(n= grid+1, model=list(ar=fi),rand.gen= function(n)sqrt(1-fi^2)*rnorm(n,0,1))
			f_grid[i,]= phi_matrix[,1:((i-1)*(grid+1))]%*%matrix(t(f_grid[(i-1):1,]),(i-1)*(grid+1),1)+ eps_grid;
		}
	}
	for (i in (p+1):tmax){
		eps_grid= arima.sim(n= grid+1, model=list(ar=fi),rand.gen= function(n)sqrt(1-fi^2)*rnorm(n,0,1))
		f_grid[i,]= phi_matrix%*%matrix(t(f_grid[(i-p):(i-1),]),p*(grid+1),1)+ eps_grid;
	}
	### The last iteration of epsilon_t is returned.
	g_cfar <- list(cfar=f_grid*sigma,epsilon=eps_grid*sigma)
}



#' Generate a CFAR(2) Process with Heteroscedasticity and Irregular Observation Locations
#'
#' Generate a convolutional functional autoregressive process of order 2 with heteroscedasticity, irregular observation locations.
#' @param tmax length of time.
#' @param grid the number of grid points used to constrcut the functional time series.
#' @param rho parameter for O-U process (noise process). 
#' @param min_obs the minimum number of observations at each time.
#' @param pois the mean for Poisson distribution. The number of observations at each follows a Poisson distribution plus min_obs.
#' @param phi_func1 the first convolutional function. Default is 0.5*x^2+0.5*x+0.13.
#' @param phi_func2 the second convolutional function. Default is 0.7*x^4-0.1*x^3-0.15*x.
#' @param weight the weight function to determine the standard deviation of O-U process (noise process). Default is 1.
#' @references
#' Liu, X., Xiao, H., and Chen, R. (2016) Convolutional autoregressive models for functional time series. \emph{Journal of Econometrics}, 194, 263-282.
#' @return The function returns a list with components:
#' \item{cfar2}{a tmax-by-(grid+1) matrix following a CFAR(1) process.}
#' \item{epsilon}{the innovation at time tmax.}
#' @examples
#' phi_func1= function(x){
#'	return(0.5*x^2+0.5*x+0.13)
#' }
#' phi_func2= function(x){
#'	return(0.7*x^4-0.1*x^3-0.15*x)
#' }
#' y=g_cfar2h(200,1000,1,40,5,phi_func1=phi_func1,phi_func2=phi_func2)
#' @export

g_cfar2h <- function(tmax=1001,grid=1000,rho=1,min_obs=40, pois=5,phi_func1=NULL, phi_func2=NULL, weight=NULL){
	#################################
	###Simulate CFAR(2) processes with heteroscedasticity, irregular observation locations
	###################################################################################

	######################################
	### We need to have f_grid to record the functional time series
	### We need to have x_pos_full to record the observation positions
	######################################

	###parameter setting
	#rho=1	### parameter for error process epsilon_t, O-U process
	#tmax= 1001;	### length of time
	#iter= 100;	### number of replications in the simulation
	#grid=1000;	### the number of grid points used to construct the functional time series X_t and epsilon_t\
	#min_obs=40	### the minimal number of observations at each time
      # phi_func1, phi_func2, weight: User specified convolution functions and weight function. The following defaults 
      #  are used if not given
	if(is.null(phi_func1)){
		phi_func1= function(x){
			return(0.5*x^2+0.5*x+0.13)
    		}
	}
	if(is.null(phi_func2)){
		phi_func2= function(x){
			return(0.7*x^4-0.1*x^3-0.15*x)
        	}
	}

	if(is.null(weight)){
		###heteroscedasticity weight function
		x_grid <- seq(0,1,by=1/grid)
		weight0= function(w){
			return((w<=0.6)*exp(-10*w)+(w>0.6)*(exp(-6)+0.2*(w-0.6))+0.1);
		}
		const=sum(weight0(x_grid)/(grid+1))
		weight= function(w){
			return(((w<=0.6)*exp(-10*w)+(w>0.6)*(exp(-6)+0.2*(w-0.6))+0.1)/const)
 		}
	}
	num_full=rpois(tmax,pois)+min_obs	###number of observations at time t follows a Poisson distribution 
                                                ###plus a number, minimal observations is required
	#########################################################################
	x_grid= seq(0, 1, by= 1/grid);	### the grid points for input variables of functional time series in [0,1]
	b_grid= seq(-1, 1, by=1/grid);	### the grid points for input variables of b-spline functions in [-1,1]
	eps_grid= matrix(0, 1, grid+1);	### the grid points for input variables of error process in [0,1]
	fi = exp(-rho/grid);		### the correlation of two adjacent points for an O-U process
	phi_b1= phi_func1(b_grid)	### phi_1(s) when s in [-1,1]
	phi_b2= phi_func2(b_grid)	### phi_2(s) when s in [-1,1]

	### the convolutional function values phi_1(s-u) and phi_2(s-u), for different s in [0,1]
	phi_matrix1= matrix(0, (grid+1), (grid+1));
	for (i in 1:(grid+1)){
		phi_matrix1[i,]= phi_b1[seq((i+grid),i,by=-1)]/(grid+1)
	}
	phi_matrix2= matrix(0, (grid+1), (grid+1));
	for (i in 1:(grid+1)){
		phi_matrix2[i,]= phi_b2[seq((i+grid),i,by=-1)]/(grid+1)	
	}

	n=max(num_full)
	x_pos_full=matrix(0,tmax,n)	###observation positions
	for(k in 1:(tmax)){
		x_pos_full[k,1:num_full[k]]=sort(runif(num_full[k]))
	}
	f_grid=matrix(0,tmax,grid+1)
	i=1
	eps_grid= arima.sim(n= grid+1, model=list(ar=fi),rand.gen= function(n)sqrt(1-fi^2)*rnorm(n,0,1))
	f_grid[1,]= eps_grid*weight(x_grid);
	i=2
	f_grid[i,]= phi_matrix1 %*% as.matrix(f_grid[i-1,])+ eps_grid*weight(x_grid);	
	for (i in 3:tmax){
		eps_grid= arima.sim(n= grid+1, model=list(ar=fi),rand.gen= function(n)sqrt(1-fi^2)*rnorm(n,0,1))
		f_grid[i,]= phi_matrix1%*%as.matrix(f_grid[i-1,]) + phi_matrix2%*%as.matrix(f_grid[i-2,])+ eps_grid*weight(x_grid);
	}

### Functional time series, number of observations at different time of periods, 
### observation points at different time of periods, the last iteration of epsilon_t is returned.
	g_cfar2h <- list(cfar2h=f_grid, num_obs=num_full, x_pos=x_pos_full,epsilon=eps_grid)
}

#' Estimation of a CFAR Process
#'
#' Estimation of a CFAR process.
#' @param f the functional time series.
#' @param p CFAR order.
#' @param df_b the degrees of freedom for natural cubic splines. Default is 10.
#' @param grid the number of gird points used to constrct the functional time series and noise process. Default is 1000.
#' @references
#' Liu, X., Xiao, H., and Chen, R. (2016) Convolutional autoregressive models for functional time series. \emph{Journal of Econometrics}, 194, 263-282.
#' @return The function returns a list with components:
#' \item{phi_coef}{estimated spline coefficients for convluaional function(s).}
#' \item{phi_func}{estimated convoluational function(s).}
#' \item{rho}{estimated rho for O-U process (noise process).}
#' \item{sigma}{estimated sigma for O-U process (noise process).}
#' @examples
#' phi_func= function(x)   
#' {
#'  	return(dnorm(x,mean=0,sd=0.1))
#' }
#' y=g_cfar1(1000,5,phi_func)
#' library(MASS)	
#' library(splines)
#' f_grid=y$cfar
#' index=seq(1,1001,by=10)
#' f=f_grid[,index]
#' est=est_cfar(f,1)
#' b_grid=seq(-1,1,by=1/grid)
#' par(mfcol=c(1,1))
#' c1 <- range(est$phi_func)
#' plot(b_grid,phi_func(b_grid),type='l',col='red',ylim=c1*1.1)
#' points(b_grid,est$phi_func,type='l')
#' @export
est_cfar <- function(f,p=3,df_b=10,grid=1000){
	if(!is.matrix(f))f <- as.matrix(f)
	t <- nrow(f)
	n=dim(f)[2]-1

	if(is.null(grid)){
		grid=1000
	}

	if(is.null(df_b)){
		df_b=10
	}

	######################################################################
	###INPUTS
	####################
	###parameter setting
	#rho=5	### parameter for error process epsilon_t, O-U process
	#t= 1000;	### length of time
	#iter= 100;	### number of replications in the simulation
	#grid=1000;	### the number of grid points used to construct the functional time series X_t and epsilon_t
	#df_b=10		### number of degrees for splines
	#n=100		### number of observations for each X_t, in the paper it is 'N'

	############################################################################
	x_grid= seq(0, 1, by= 1/grid);	### the grid points for input variables of functional time series in [0,1]
	b_grid= seq(-1, 1, by=1/grid);	### the grid points for input variables of b-spline functions in [-1,1]
	eps_grid= matrix(0, 1, grid+1);	### the grid points for input variables of error process in [0,1]
	#############################
	### It shows the best approximation of phi using b-spline with df=df_b

	coef_b=df_b+1;
	b_grid_sp= ns(b_grid, df=df_b);

	###############################
	###Estimation
	x= seq(0, 1, by=1/n);
	index= 1:(n+1)*(grid/n)-(grid/n)+1
	x_grid_sp= bs(x_grid, df=n, degree=1);	###basis function values if we interpolate discrete observations of x_t 
	x_sp= x_grid_sp[index,];			###basis function values at each observation point 
	df=n;
	coef= df+1

	###convolutions of b-pline functions (df=df) for phi_1 and phi_2 and basis functions (df=n) for x
	x_grid_full= cbind(rep(1,grid+1), x_grid_sp);
	b_grid_full= cbind(rep(1,2*grid+1), b_grid_sp);
	b_matrix=matrix(0,grid+1,grid+1)
	bstar_grid2= matrix(0, grid+1, coef*coef_b)
	for (j in 1:(coef_b)){
		for (i in 1:(grid+1)){
			b_matrix[i,]= b_grid_full[seq((i+grid),i, by=-1),j]/(grid+1);
		}
		bstar_grid2[, (j-1)*coef+(1:coef)]= b_matrix %*% x_grid_full
	}

	bstar2= bstar_grid2[index,]
	bstar3= matrix(0, (n+1)*coef_b, coef)
	for (i in 1:coef_b){
		bstar3[(i-1)*(n+1)+(1:(n+1)),]=bstar2[, (i-1)*coef+(1:coef)]
	}


	indep=matrix(0,(n+1)*(t-1),coef_b)	### design matrix for b-spline coefficient of phi_1
	dsg=cbind(rep(1,n+1),x_sp);
	dsg_mat= solve(t(dsg)%*%dsg)%*%t(dsg)
	ahat= f%*%t(dsg_mat)	### b-spline coefficient when we interpolate x_t
	for (i in 2:t){	
		M_t= matrix(ahat[i-1,]%*% t(bstar3), nrow=n+1, ncol=df_b+1)
		indep[(i-2)*(n+1)+(1:(n+1)),]=M_t	### design matrix for b-spline coefficient of phi_1
	}

	indep.tmp=NULL
	for(i in 1:p){
		tmp=indep[((i-1)*(n+1)+1):((n+1)*(t-p-1+i)),]
		indep.tmp=cbind(tmp,indep.tmp)
	}
	indep.p=indep.tmp

	pdf4=function(para4)
	{	
		psi= para4;	### correlation between two adjacent observation point exp(-rho/n)
		psi_matinv= matrix(0, n+1, n+1)	### correlation matrix of error process at observation points
		psi_matinv[1,1]=1
		psi_matinv[n+1,n+1]=1;
		psi_matinv[1,2]=-psi;
		psi_matinv[n+1,n]=-psi;
		for (i in 2:n){
			psi_matinv[i,i]= 1+psi^2;
			psi_matinv[i,i+1]= -psi;
			psi_matinv[i,i-1]= -psi;
		}
		mat1= matrix(0, coef_b*p, coef_b*p);
		mat2= matrix(0, coef_b*p, 1)
		for(i in (p+1):t){
			tmp.n=n+1
			M_t=indep.p[(i-p-1)*(n+1)+(1:tmp.n),]
			tmp_mat= t(M_t) %*% psi_matinv;
			mat1= mat1+ tmp_mat %*% M_t;
			mat2= mat2+ tmp_mat %*% f[i,];
		}
		phihat= solve(mat1)%*%mat2
		ehat=0
		log.mat=0
		for (i in (p+1):t){
			tmp.n=n+1
			M_t= indep.p[(i-p-1)*(n+1)+(1:tmp.n),]
			eps= f[i,1:tmp.n]- M_t%*% phihat
			epart= t(eps)%*%psi_matinv %*%eps
			ehat= epart+ehat
			log.mat=log.mat+ log(det(psi_matinv))
		}
		l=-(t-p)*n/2*log(ehat)+1/2*log.mat
		return(-l);
	}
	para4= exp(-5)
	result4=optim(para4,pdf4, lower=0.001, upper=0.9999, method='L-BFGS-B')
	psi=result4$par
	psi_matinv= matrix(0, n+1, n+1)
	psi_matinv[1,1]=1
	psi_matinv[n+1,n+1]=1;
	psi_matinv[1,2]=-psi;
	psi_matinv[n+1,n]=-psi;
	for (i in 2:n){
		psi_matinv[i,i]= 1+psi^2;
		psi_matinv[i,i+1]= -psi;
		psi_matinv[i,i-1]= -psi;
	}
	mat1= matrix(0, p*(df_b+1), p*(df_b+1));	### matrix under the first brackets in formula (15)
	mat2= matrix(0, p*(df_b+1), 1);	### matrix under the second brackets in formula (15)
	for (i in (p+1):t){	
		M_t= indep.p[(i-p-1)*(n+1)+(1:(n+1)),]
		tmp_mat= t(M_t) %*% psi_matinv;
		mat1= mat1+ tmp_mat %*% M_t;
		mat2= mat2+ tmp_mat %*% f[i,];
	}
	phihat= solve(mat1)%*%mat2
	phihat_func=matrix(0,p,(1+2*grid))
	for(i in 1:p){
		phihat_func[i,]= cbind(rep(1,2*grid+1),b_grid_sp)%*%phihat[(1+(coef_b*(i-1))):(i*coef_b)];		###b-spline coefficient of estimated phi_1
	}
	ehat=0
	for (i in (p+1):t){
		tmp.n= n+1
		M_t= indep.p[(i-p-1)*(n+1)+(1:tmp.n),]
		eps= f[i,1:tmp.n]- M_t%*% phihat
		epart= t(eps)%*%psi_matinv %*%eps
		ehat= epart+ehat
	}
	rho_hat= -log(psi)*n
	sigma_hat=ehat/((t-p)*(n+1))*2*rho_hat/(1-psi^2)
	phi_coef=t(matrix(phihat,coef_b,p))
	est_cfar <- list(phi_coef=phi_coef,phi_func=phihat_func,rho=rho_hat,sigma2=sigma_hat)
	return(est_cfar)
}


#' F Test for a CFAR Process
#'
#' F test for a CFAR process to specify CFAR order.
#' @param f the functional time series.
#' @param p.max maximum CFAR order. Default is 6.
#' @param df_b the degrees of freedom for natural cubic splines. Default is 10.
#' @param grid the number of gird points used to constrct the functional time series and noise process. Default is 1000.
#' @references
#' Liu, X., Xiao, H., and Chen, R. (2016) Convolutional autoregressive models for functional time series. \emph{Journal of Econometrics}, 194, 263-282.
#' @return The function outputs F test statistics and their p-values.
#' @examples
#' phi_func= function(x)   
#' {
#'  	return(dnorm(x,mean=0,sd=0.1))
#' }
#' y=g_cfar1(1000,5,phi_func)
#' library(MASS)	
#' library(splines)
#' f_grid=y$cfar
#' index=seq(1,1001,by=10)
#' f=f_grid[,index]
#' F_test_cfar(f)
#' @export
F_test_cfar <- function(f,p.max=6,df_b=10,grid=1000){
#	library(MASS)	
#	library(splines)
	##########################
	### F test: CFAR(0)vs CFAR(1), CFAR(1)vs CFAR(2), CFAR(2)vs CFAR(3) when rho is known.
	#
	if(!is.matrix(f))f <- as.matrix(f)
	t <- nrow(f)
	n=dim(f)[2]-1

	if(is.null(p.max)){
		p.max=6
	}
	if(is.null(df_b)){
		df_b=10
	}
	if(is.null(grid)){
		grid=1000
	}
	##################################################################################
	###parameter setting
	#rho=5	### parameter for error process epsilon_t, O-U process
	#t= 1000;	### length of time
	##iter= 100;	### number of replications in the simulation
	#grid=1000;	### the number of grid points used to construct the functional time series X_t and epsilon_t
	#df_b=10		### number of degrees for splines
	#n=100		### number of observations for each X_t, in the paper it is 'N'
	###############################################################################
	x_grid= seq(0, 1, by= 1/grid);	### the grid points for input variables of functional time series in [0,1]
	b_grid= seq(-1, 1, by=1/grid);	### the grid points for input variables of b-spline functions in [-1,1]
	eps_grid= matrix(0, 1, grid+1);	### the grid points for input variables of error process in [0,1]
	coef_b=df_b+1;				### number of df for b-spline
	b_grid_sp= ns(b_grid, df=df_b);	###b-spline functions
	x= seq(0, 1, by=1/n);			###observations points for X_t process
	index= 1:(n+1)*(grid/n)-(grid/n)+1
	x_grid_sp= bs(x_grid, df=n, degree=1);	###basis function values if we interpolate discrete observations of x_t
	x_sp= x_grid_sp[index,];			###basis function values at each observation point
	df=n;			### number of interplolation of X_t
	coef= df+1;
	###convolutions of b-pline functions (df=df_b) for phi and basis functions (df=n) for x
	x_grid_full= cbind(rep(1,grid+1), x_grid_sp);
	b_grid_full= cbind(rep(1,2*grid+1), b_grid_sp);
	b_matrix=matrix(0,grid+1,grid+1)
	bstar_grid2= matrix(0, grid+1, coef*coef_b)
	for (j in 1:(coef_b)){
		for (i in 1:(grid+1)){
			b_matrix[i,]= b_grid_full[seq((i+grid),i, by=-1),j]/(grid+1);
		}
		bstar_grid2[, (j-1)*coef+(1:coef)]= b_matrix %*% x_grid_full
	###	print(j)
	}
	bstar2= bstar_grid2[index,]
	bstar3= matrix(0, (n+1)*coef_b, coef)
	for (i in 1:coef_b){
		bstar3[(i-1)*(n+1)+(1:(n+1)),]=bstar2[, (i-1)*coef+(1:coef)]
	}
	dsg=cbind(rep(1,n+1),x_sp);
	dsg_mat= solve(t(dsg)%*%dsg)%*%t(dsg)
	ahat= f%*%t(dsg_mat)
	indep=matrix(0,(n+1)*(t-1),coef_b)
	for (i in 2:t){	
		M_t= matrix(ahat[i-1,]%*% t(bstar3), nrow=n+1, ncol=df_b+1)
		indep[(i-2)*(n+1)+(1:(n+1)),]=M_t
	}
#########################################
	######### AR(0) and AR(1)
	#########################################
	#######################
	####cat("Data set: ",q,"\n")
	### Fit cfar(1) model
	###	f= f_grid[(q-1)*tmax+(1:t)+(tmax-t), index];
	p=1
	dsg=cbind(rep(1,n+1),x_sp);
	dsg_mat= solve(t(dsg)%*%dsg)%*%t(dsg)
	ahat= f%*%t(dsg_mat)
	indep=matrix(0,(n+1)*(t-1),coef_b)
	for (i in 2:t){	
		M_t= matrix(ahat[i-1,]%*% t(bstar3), nrow=n+1, ncol=df_b+1)
		indep[(i-2)*(n+1)+(1:(n+1)),]=M_t
	}
	est=est_cfar(f,1,df_b=df_b,grid)
	phihat= est$phi_coef
	psi=exp(-est$rho/n)
	mat1= matrix(0, coef_b*p, coef_b*p);
	mat2= matrix(0, coef_b*p, 1);
	phi_matinv=matrix(0, n+1, n+1)
	phi_matinv[1,1]=1;
	phi_matinv[n+1,n+1]=1;
	phi_matinv[1,2]= -psi;
	phi_matinv[n+1,n]= -psi
	for (i in 2:n){
		phi_matinv[i,i]= 1+psi^2;
		phi_matinv[i,i+1]= -psi;
		phi_matinv[i,i-1]= -psi
	}

	predict= t(matrix(indep%*%matrix(phihat,df_b+1,1),ncol=t-1,nrow=n+1))
	ssr1= sum(diag((predict)%*%phi_matinv%*% t(predict)))
	sse1= sum(diag((f[2:t,]-predict[1:(t-1),])%*%phi_matinv%*% t(f[2:t,]-predict[1:(t-1),])))
	sse0=sum(diag(f[2:t,]%*%phi_matinv%*% t(f[2:t,])))
	statistic= (sse0-sse1)/coef_b/sse1*((n+1)*(t-2)-coef_b)
	pval=1-pf(statistic,coef_b,(t-2)*(n+1)-coef_b)
	cat("Test and p-value of Order 0 vs Order 1: ","\n")
	print(c(statistic,pval))

	p=1
	sse.pre= sum(diag((f[(p+2):t,]-predict[2:(t-p),])%*%phi_matinv%*% t(f[(p+2):t,]-predict[2:(t-p),])))

	for(p in 2:p.max){
		indep.tmp=NULL
		for(i in 1:p){
			tmp=indep[((i-1)*(n+1)+1):((n+1)*(t-p-1+i)),]
			indep.tmp=cbind(tmp,indep.tmp)
		}
		indep.p=indep.tmp
		pdf4=function(para4)
		{	
			mat1= matrix(0, coef_b*p, coef_b*p);
			mat2= matrix(0, coef_b*p, 1);
			psi= para4;
			phi_matinv=matrix(0, n+1, n+1)
			phi_matinv[1,1]=1;
			phi_matinv[n+1,n+1]=1;
			phi_matinv[1,2]= -psi;
			phi_matinv[n+1,n]= -psi
			for (i in 2:n){
				phi_matinv[i,i]= 1+psi^2;
				phi_matinv[i,i+1]= -psi;
				phi_matinv[i,i-1]= -psi
			}
			psi_matinv=phi_matinv
			for(i in (p+1):t){
				tmp.n=n+1
				tmp.index=which(!is.na(f[i,]))
				M_t=indep.p[(i-p-1)*(n+1)+(1:tmp.n),]
				tmp_mat= t(M_t) %*% psi_matinv;
				mat1= mat1+ tmp_mat %*% M_t;
				mat2= mat2+ tmp_mat %*% f[i,1:tmp.n];
			}
			phihat= solve(mat1)%*%mat2
			ehat=0
			log.mat=0
			for (i in (p+1):t){
				tmp.n=n+1
				M_t= indep.p[(i-p-1)*(n+1)+(1:tmp.n),]
				eps= f[i,1:tmp.n]- M_t%*% phihat
				epart= t(eps)%*%psi_matinv %*%eps
				ehat= epart+ehat
				log.mat=log.mat+ log(det(psi_matinv))
			}
			l=-((t-p)*(n+1))/2*log(ehat)+1/2*log.mat
			return(-l);
		}
		para4= exp(-5)
		result4=optim(para4,pdf4, lower=0.001, upper=0.999, method='L-BFGS-B')

		mat1= matrix(0, p*coef_b, p*coef_b);
		mat2= matrix(0, coef_b*p, 1);
		psi_invall=matrix(0,(n+1)*(t-1),(n+1))
		psi=result4$par
		phi_matinv=matrix(0, n+1, n+1)
		phi_matinv[1,1]=1;
		phi_matinv[n+1,n+1]=1;
		phi_matinv[1,2]= -psi;
		phi_matinv[n+1,n]= -psi
		psi_matinv=phi_matinv
		for (i in 2:n){
			phi_matinv[i,i]= 1+psi^2;
			phi_matinv[i,i+1]= -psi;
			phi_matinv[i,i-1]= -psi
		}
		psi_matinv=phi_matinv
		for(i in (p+1):t){
			tmp.n=n+1
			M_t= indep.p[(i-p-1)*(n+1)+(1:tmp.n),]
			tmp_mat= t(M_t) %*% psi_matinv;
			mat1= mat1+ tmp_mat %*% M_t;
			mat2= mat2+ tmp_mat %*% f[i,1:tmp.n];
		}
		phihat= solve(mat1)%*%mat2
		predict= t(matrix(indep.p%*%phihat,ncol=t-p,nrow=n+1))
		ssr.p= sum(diag((predict)%*%phi_matinv%*% t(predict)))
		sse.p= sum(diag((f[(p+1):t,]-predict)%*%phi_matinv%*% t(f[(p+1):t,]-predict)))
		statistic= (sse.pre-sse.p)/coef_b/sse.p*((n+1)*(t-p-1)-p*coef_b)
		pval=1-pf(statistic,coef_b,(t-p)*(n+1)-p*coef_b)
		cat("Test and  p-value of Order", p-1, "vs Order",p,": ","\n")
		print(c(statistic,pval))
		sse.pre= sum(diag((f[(p+2):t,]-predict[2:(t-p),])%*%phi_matinv%*% t(f[(p+2):t,]-predict[2:(t-p),])))
	}
}



#' F Test for a CFAR Process with Heteroscedasticity and Irregular Observation Locations
#'
#' F test for a CFAR process with heteroscedasticity and irregular observation locations to specify the CFAR order.
#' @param f the functional time series.
#' @param weight the covariance functions for noise process.
#' @param p.max maximum CFAR order. Default is 3.
#' @param grid the number of gird points used to constrct the functional time series and noise process. Default is 1000.
#' @param df_b the degrees of freedom for natural cubic splines. Default is 10.
#' @param num_obs the numbers of observations. It is a t-by-1 vector, where t is the length of time.
#' @param x_pos the observation location matrix. If the locations are regular, it is a t-by-(n+1) matrix with all entries 1/n.
#' @references
#' Liu, X., Xiao, H., and Chen, R. (2016) Convolutional autoregressive models for functional time series. \emph{Journal of Econometrics}, 194, 263-282.
#' @return The function outputs F test statistics and their p-values.
#' @examples
#' phi_func1= function(x){
#'	return(0.5*x^2+0.5*x+0.13)
#' }
#' phi_func2= function(x){
#'	return(0.7*x^4-0.1*x^3-0.15*x)
#' }
#' library(MASS)	
#' library(splines)
#' testh=g_cfar2h(100,1000,1,40,5,phi_func1=phi_func1,phi_func2=phi_func2)
#' #obtain discret observations
#' f_grid=testh$cfar2h
#' num_obs=testh$num_obs
#' x_pos=testh$x_pos
#' t=dim(f_grid)[1]
#' n=max(num_obs)
#' index=matrix(0,t,n)
#' x_grid=seq(0,1,by=1/grid)
#' f=matrix(0,t,n)
#' for(i in 1:t){
#'	for(j in 1:num_obs[i]){
#'		index[i,j]=which(abs(x_pos[i,j]-x_grid)==min(abs(x_pos[i,j]-x_grid)))
#'	}
#'	f[i,1:num_obs[i]]=f_grid[i,index[i,1:num_obs[i]]];
#' }
#' x_grid <- seq(0,1,by=1/grid)
#' weight0= function(w){
#'	return((w<=0.6)*exp(-10*w)+(w>0.6)*(exp(-6)+0.2*(w-0.6))+0.1);
#' }
#' const=sum(weight0(x_grid)/(grid+1))
#' weight= function(w){
#'	return(((w<=0.6)*exp(-10*w)+(w>0.6)*(exp(-6)+0.2*(w-0.6))+0.1)/const)
#' }
#' F_test_cfarh(f,weight,3,1000,5,num_obs,x_pos)
#' @export
F_test_cfarh <- function(f,weight,p.max=3,grid=1000,df_b=10,num_obs=NULL,x_pos=NULL){
#	library(MASS)	
#	library(splines)
	########################

	if(!is.matrix(f))f <- as.matrix(f)
	t <- nrow(f)

	if(is.null(p.max)){
		p.max=3
	}
	
	if(is.null(num_obs)){
		num_obs=dim(f)[2]
		n=dim(f)[2]
	}else{
		n=max(num_obs)
	}

	if(length(num_obs)!=t){
		num_obs=rep(n,t)
	}
	if(is.null(df_b)){
		df_b=10
	}
	if(is.null(grid)){
		grid=1000
	}
	if(is.null(x_pos)){
		x_pos=matrix(rep(seq(0,1,by=1/(num_obs[1]-1)),each=t),t,num_obs[1])
	}

	##################
	###parameter setting
	#grid=1000;	### the number of grid points used to construct the functional time series X_t and epsilon_t\
	#coef_b=6	### k=6
	# num_obs is a t by 1 vector which records N_t for each time

	#######################################################K
	x_grid=seq(0,1,by=1/grid)
	coef_b=df_b+1
	b_grid=seq(-1,1,by=1/grid)
	b_grid_sp = ns(b_grid,df=df_b);
	b_grid_full= cbind(rep(1,2*grid+1), b_grid_sp);
	indep= matrix(0, (n+1)*(t-1),coef_b)
	bstar_grid= matrix(0, grid+1, coef_b)
	b_matrix= matrix(0, (grid+1), (grid+1))
	bstar_grid_full= matrix(0, (grid+1)*t,coef_b)
	index=matrix(0,t,n)
	for(i in 1:t){
		for(j in 1:num_obs[i]){
			index[i,j]=which(abs(x_pos[i,j]-x_grid)==min(abs(x_pos[i,j]-x_grid)))
		}
	}
	for (i in 1:(t-1)){
		rec_x= approx(x_pos[i,1:num_obs[i]],f[i,1:num_obs[i]],xout=x_grid,rule=2,method='linear')	
                         ### rec_x is the interpolation of x_t
		for(j in 1:coef_b){
			for (k in 1:(grid+1)){
				b_matrix[k,]= b_grid_full[seq((k+grid),k, by=-1),j]/(grid+1);
			}
			bstar_grid[, j]= b_matrix %*% matrix(rec_x$y,ncol=1,nrow=grid+1)
		}
		bstar_grid_full[(i-1)*(grid+1)+1:(grid+1),]=bstar_grid	###convoluttion of basis spline function and x_t
		tmp=bstar_grid[index[i+1,1:num_obs[i+1]],]
		indep[(i-1)*(n+1)+1:num_obs[i+1],]=tmp
	}



	
	#################AR(1)
	pdf4=function(para4)	###Q function in formula (12)
	{	
		mat1= matrix(0, df_b+1, df_b+1);
		mat2= matrix(0, df_b+1, 1);
		### correlation matrix of error process at observation points
		psi_invall=matrix(0,(n+1)*(t-1),(n+1))
		for(i in 2:t){
			psi= para4;
			tmp.n=num_obs[i]
			psi_mat= matrix(1, tmp.n, tmp.n)
			for(k in 1:(tmp.n)){
				for(j in 1:(tmp.n)){
					psi_mat[k,j]= psi^(abs(x_pos[i,j]-x_pos[i,k]))
				}
			}
			psi_matinv2=solve(diag(weight(x_pos[i,1:tmp.n]))%*%psi_mat%*%diag(weight(x_pos[i,1:tmp.n])))
			psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]=psi_matinv2
			M_t=indep[(i-2)*(n+1)+(1:tmp.n),]
			tmp_mat= t(M_t) %*% psi_matinv2;
			mat1= mat1+ tmp_mat %*% M_t;			### matrix under the first brackets in formula (15)
			mat2= mat2+ tmp_mat %*% f[i,1:tmp.n];	### matrix under the second brackets in formula (15)
		}
		phihat= solve(mat1)%*%mat2
		ehat=0	###e(beta,phi) in formula (12)
		log.mat=0
		for (i in 2:t){
			tmp.n=num_obs[i]
			M_t= indep[(i-2)*(n+1)+(1:tmp.n),]
			eps= f[i,1:tmp.n]- M_t%*% phihat
			psi_matinv= psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]
			epart= t(eps)%*%psi_matinv %*%eps
			ehat= epart+ehat
			log.mat=log.mat+ log(det(psi_matinv))
		}
		l=-sum(num_obs[2:t])/2*log(ehat)+1/2*log.mat	### the number defined in formula (14)
		return(-l);
	}
	para4= exp(-1)
	result4=optim(para4,pdf4, lower=0.001, upper=0.999, method='L-BFGS-B')
	mat1= matrix(0, df_b+1, df_b+1);
	mat2= matrix(0, df_b+1, 1);
	psi_invall=matrix(0,(n+1)*(t-1),(n+1))
	psi=result4$par
	for(i in 2:t){
		tmp.n=num_obs[i]
		psi_mat= matrix(1, tmp.n, tmp.n)
		for(k in 1:tmp.n){
			for(j in 1:tmp.n){
				psi_mat[k,j]=psi^(abs(x_pos[i,j]-x_pos[i,k]))
			}
		}
		psi_matinv2=solve(diag(weight(x_pos[i,1:tmp.n]))%*%psi_mat%*%diag(weight(x_pos[i,1:tmp.n])))	
		psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]= psi_matinv2
		M_t= indep[(i-2)*(n+1)+(1:tmp.n),]
		tmp_mat= t(M_t) %*% psi_matinv2;
		mat1= mat1+ tmp_mat %*% M_t;
		mat2= mat2+ tmp_mat %*% f[i,1:tmp.n];
	}
	phihat= solve(mat1)%*%mat2	
	ehat=0		###e(beta,phi) in formula (12)
	for (i in 2:t){
		tmp.n=num_obs[i]
		M_t= indep[(i-2)*(n+1)+(1:tmp.n),]
		eps= f[i,1:tmp.n]- M_t%*% phihat
		psi_matinv=psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]
		epart= t(eps)%*%psi_matinv %*%eps
		ehat= epart+ehat
	}
	rho_hat= -log(psi)
	sigma_hat=ehat/sum(num_obs[2:t])*2*rho_hat

	sigma_hat1=sigma_hat
	rho_hat1=rho_hat
	predict1= t(matrix(indep%*%phihat,ncol=t-1,nrow=n+1))
	ssr1=0
	sse1=0
	for (i in 2:t){
		tmp.n=num_obs[i]
		psi_matinv= psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]
		ssr1= ssr1+ sum((predict1[i-1,1:tmp.n] %*% psi_matinv) * (predict1[i-1,1:tmp.n]))
		sse1= sse1+ sum(((f[i,1:tmp.n]- predict1[i-1,1:tmp.n]) %*% psi_matinv) *(f[i,1:tmp.n]-predict1[i-1,1:tmp.n]))
	}

###############################
	######   No AR term
	psi=exp(-rho_hat1)
	mat1= matrix(0, df_b+1, df_b+1);
	mat2= matrix(0, df_b+1, 1);
	phihat=matrix(0,coef_b,1)
	ehat=0
	for (i in 2:t){
		tmp.n=num_obs[i]
		M_t= indep[(i-2)*(n+1)+(1:tmp.n),]
		eps= f[i,1:tmp.n]- M_t%*% phihat
		psi_matinv=psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]
		epart= t(eps)%*%psi_matinv %*%eps
		ehat= epart+ehat
	}
	rho_hat= -log(psi)
	sigma_hat=ehat/sum(num_obs[2:t])*2*rho_hat
	rho_hat0=rho_hat
	sigma_hat0=sigma_hat
	sse0=0
	test=0
	for (i in 2:t){
		tmp.n=num_obs[i]
		psi_matinv= psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]
		sse0= sse0+ sum((f[i,1:tmp.n] %*% psi_matinv) *(f[i,1:tmp.n]))
	}

	statistic= (sse0-sse1)/coef_b/sse1*(sum(num_obs[2:t])-coef_b)
	pval=1-pf(statistic,coef_b,sum(num_obs[2:t])-coef_b)
	cat("Test and  p-value of Order 0 vs Order 1: ","\n")
	print(c(statistic,pval))
	
	
	indep.p=indep
	if (p.max>1){
		for(p in 2:p.max){
			indep.pre=indep.p[(n+2):((n+1)*(t-p+1)),]
			indep.tmp=indep
			for (i in 1:(t-p)){
				for(j in 1:num_obs[i+p]){
					index[i+p,j]= which(abs(x_pos[i+p,j]-x_grid)==min(abs(x_pos[i+p,j]-x_grid)))
				}
				tmp=bstar_grid_full[(i-1)*(grid+1)+index[i+p,1:num_obs[i+p]],]
				indep.tmp[(i-1)*(n+1)+1:num_obs[i+p],]=tmp
			}
			indep.test=cbind(indep.p[(n+2):((n+1)*(t-p+1)),], indep.tmp[1:((n+1)*(t-p)),])
			indep.p=indep.test

			pdf4=function(para4)	###Q function in formula (12)
			{	
				mat1= matrix(0, p*(df_b+1), p*(df_b+1));
				mat2= matrix(0, p*(df_b+1), 1);
				### correlation matrix of error process at observation points
				psi_invall=matrix(0,(n+1)*(t-1),(n+1))
				for(i in (p+1):t){
					psi= para4;
					tmp.n=num_obs[i]
					psi_mat= matrix(1, tmp.n, tmp.n)
					for(k in 1:(tmp.n)){
						for(j in 1:(tmp.n)){
							psi_mat[k,j]= psi^(abs(x_pos[i,j]-x_pos[i,k]))
						}
					}
					psi_matinv2=solve(diag(weight(x_pos[i,1:tmp.n]))%*%psi_mat%*%diag(weight(x_pos[i,1:tmp.n])))
					psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]=psi_matinv2
					M_t=indep.p[(i-p-1)*(n+1)+(1:tmp.n),]
					tmp_mat= t(M_t) %*% psi_matinv2;
					mat1= mat1+ tmp_mat %*% M_t;			### matrix under the first brackets in formula (15)
					mat2= mat2+ tmp_mat %*% f[i,1:tmp.n];	### matrix under the second brackets in formula (15)
				}
				phihat= solve(mat1)%*%mat2
				ehat=0	###e(beta,phi) in formula (12)
				log.mat=0
				for (i in (p+1):t){
					tmp.n=num_obs[i]
					M_t= indep.p[(i-p-1)*(n+1)+(1:tmp.n),]	
					eps= f[i,1:tmp.n]- M_t%*% phihat
					psi_matinv= psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]
					epart= t(eps)%*%psi_matinv %*%eps
					ehat= epart+ehat
					log.mat=log.mat+ log(det(psi_matinv))
				}
				l=-sum(num_obs[(p+1):t])/2*log(ehat)+1/2*log.mat	### the number defined in formula (14)
				return(-l);
			}
			para4= exp(-1)
			result4=optim(para4,pdf4, lower=0.001, upper=0.999, method='L-BFGS-B')
			mat1= matrix(0, p*(df_b+1), p*(df_b+1));
			mat2= matrix(0, p*(df_b+1), 1);
			psi_invall=matrix(0,(n+1)*(t-1),(n+1))
			psi=result4$par
			for(i in (p+1):t){
				tmp.n=num_obs[i]
				psi_mat= matrix(1, tmp.n, tmp.n)
				for(k in 1:tmp.n){
					for(j in 1:tmp.n){
						psi_mat[k,j]=psi^(abs(x_pos[i,j]-x_pos[i,k]))
					}
				}
				psi_matinv2=solve(diag(weight(x_pos[i,1:tmp.n]))%*%psi_mat%*%diag(weight(x_pos[i,1:tmp.n])))	
				psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]= psi_matinv2
				M_t= indep.p[(i-p-1)*(n+1)+(1:tmp.n),]
				tmp_mat= t(M_t) %*% psi_matinv2;
				mat1= mat1+ tmp_mat %*% M_t;
				mat2= mat2+ tmp_mat %*% f[i,1:tmp.n];
			}
			phihat= solve(mat1)%*%mat2	
			ehat=0		###e(beta,phi) in formula (12)
			for (i in (p+1):t){
				tmp.n=num_obs[i]
				M_t= indep.p[(i-p-1)*(n+1)+(1:tmp.n),]
				eps= f[i,1:tmp.n]- M_t%*% phihat
				psi_matinv=psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]
				epart= t(eps)%*%psi_matinv %*%eps
				ehat= epart+ehat
			}
			rho_hat= -log(psi)
			sigma_hat=ehat/sum(num_obs[(p+1):t])*2*rho_hat

			sigma_hat1=sigma_hat
			rho_hat1=rho_hat
			predict= t(matrix(indep.p%*%phihat,ncol=t-p,nrow=n+1))
			ssr.p= 0
			sse.p=0
			for (i in (p+1):t){
				tmp.n=num_obs[i]
				psi_matinv= psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]
				ssr.p= ssr.p+ sum((predict[i-p,1:tmp.n] %*% psi_matinv) * (predict[i-p,1:tmp.n]))
				sse.p= sse.p+ sum(((f[i,1:tmp.n]- predict[i-p,1:tmp.n]) %*% psi_matinv) *(f[i,1:tmp.n]-predict[i-p,1:tmp.n]))
			}

			mat1= matrix(0, coef_b*(p-1), coef_b*(p-1));
			mat2= matrix(0, coef_b*(p-1), 1);
			psi_invall=matrix(0,(n+1)*(t-1),(n+1))
			for(i in (p+1):t){
				tmp.n=num_obs[i]
				tmp.index=which(!is.na(f[i,]))
				psi_mat= matrix(1, tmp.n, tmp.n)
				for(k in 1:tmp.n){
					for(j in 1:tmp.n){
						psi_mat[k,j]=psi^(abs(x_pos[i,j]-x_pos[i,k]))
					}
				}
				psi_matinv2=solve(diag(weight(x_pos[i,1:tmp.n]))%*%psi_mat%*%diag(weight(x_pos[i,1:tmp.n])))
				psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]= psi_matinv2
				M_t= indep.pre[(i-p-1)*(n+1)+(1:tmp.n),]
				tmp_mat= t(M_t) %*% psi_matinv2;
				mat1= mat1+ tmp_mat %*% M_t;
				mat2= mat2+ tmp_mat %*% f[i,1:tmp.n];
			}
			phihat= solve(mat1)%*%mat2
			predict= t(matrix(indep.pre%*%phihat,ncol=t-p,nrow=n+1))
			ssr.pre=0
			sse.pre=0.0000
			for (i in (p+1):t){
				tmp.n=num_obs[i]
				tmp.index=which(!is.na(f[i,]))
				psi_matinv= psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]
				ssr.pre= ssr.pre+ sum((predict[i-p,1:tmp.n] %*% psi_matinv) * (predict[i-p,1:tmp.n]))
				sse.pre= sse.pre+ sum(((f[i,1:tmp.n]- predict[i-p,1:tmp.n]) %*% psi_matinv) *(f[i,1:tmp.n]-predict[i-p,1:tmp.n]))
			}
			statistic= (sse.pre-sse.p)/coef_b/sse.p*(sum(num_obs[(p+1):t])-coef_b)
			pval=1-pf(statistic,coef_b,sum(num_obs[(p+1):t])-coef_b)
			cat("Test and  p-value of Order",p-1," vs Order", p,": ","\n")
			print(c(statistic,pval))	
		}
	}
}




#' Estimation of a CFAR Process with Heteroscedasticity and Irregualr Observation Locations
#'
#' Estimation of a CFAR process with heteroscedasticity and irregualr observation locations.
#' @param f the functional time series.
#' @param weight the covariance functions of noise process.
#' @param p CFAR order.
#' @param grid the number of gird points used to constrct the functional time series and noise process. Default is 1000.
#' @param df_b the degrees of freedom for natural cubic splines. Default is 10.
#' @param num_obs the numbers of observations. It is a t-by-1 vector, where t is the length of time.
#' @param x_pos the observation location matrix. If the locations are regular, it is a t-by-(n+1) matrix with all entries 1/n.
#' @references
#' Liu, X., Xiao, H., and Chen, R. (2016) Convolutional autoregressive models for functional time series. \emph{Journal of Econometrics}, 194, 263-282.
#' @return The function returns a list with components:
#' \item{phi_coef}{estimated spline coefficients for convluaional function(s).}
#' \item{phi_func}{estimated convoluational function(s).}
#' \item{rho}{estimated rho for O-U process (noise process).}
#' \item{sigma}{estimated sigma for O-U process (noise process).}
#' @examples
#' phi_func1= function(x){
#'	return(0.5*x^2+0.5*x+0.13)
#' }
#' phi_func2= function(x){
#'	return(0.7*x^4-0.1*x^3-0.15*x)
#' }
#' library(MASS)	
#' library(splines)
#' testh=g_cfar2h(100,1000,1,40,5,phi_func1=phi_func1,phi_func2=phi_func2)
#' #obtain discret observations
#' f_grid=testh$cfar2h
#' num_obs=testh$num_obs
#' x_pos=testh$x_pos
#' t=dim(f_grid)[1]
#' n=max(num_obs)
#' index=matrix(0,t,n)
#' x_grid=seq(0,1,by=1/grid)
#' f=matrix(0,t,n)
#' for(i in 1:t){
#'	for(j in 1:num_obs[i]){
#'		index[i,j]=which(abs(x_pos[i,j]-x_grid)==min(abs(x_pos[i,j]-x_grid)))
#'	}
#'	f[i,1:num_obs[i]]=f_grid[i,index[i,1:num_obs[i]]];
#' }
#' x_grid <- seq(0,1,by=1/grid)
#' weight0= function(w){
#'	return((w<=0.6)*exp(-10*w)+(w>0.6)*(exp(-6)+0.2*(w-0.6))+0.1);
#' }
#' const=sum(weight0(x_grid)/(grid+1))
#' weight= function(w){
#'	return(((w<=0.6)*exp(-10*w)+(w>0.6)*(exp(-6)+0.2*(w-0.6))+0.1)/const)
#' }
#' est=est_cfarh(f,weight,2,1000,5,num_obs,x_pos)
#' par(mfrow=c(1,2))
#' c1 <- range(est$phi_func)
#' plot(b_grid,phi_func1(b_grid),type='l',col='red',ylim=c1*1.05)
#' points(b_grid,est$phi_func[1,],type='l')
#' c1 <- range(est$phi_func)
#' plot(b_grid,phi_func2(b_grid),type='l',col='red',ylim=c1*1.05)
#' points(b_grid,est$phi_func[2,],type='l')
#' cat("rho: ",est$rho,"\n")
#' cat("sigma2: ",est$sigma2,"\n")
#' @export

est_cfarh <- function(f,weight,p=2,grid=1000,df_b=5, num_obs=NULL,x_pos=NULL){
#	library(MASS)	
#	library(splines)
	########################
	### CFAR(2) with processes with heteroscedasticity, irregular observation locations
	### Estimation of phi(), rho and sigma, and F test 

	if(!is.matrix(f))f <- as.matrix(f)
	t <- nrow(f)
	
	if(is.null(num_obs)){
		num_obs=dim(f)[2]
		n=dim(f)[2]
	}else{
		n=max(num_obs)
	}

	if(length(num_obs)!=t){
		num_obs=rep(n,t)
	}

	if(is.null(x_pos)){
		x_pos=matrix(rep(seq(0,1,by=1/(num_obs[1]-1)),each=t),t,num_obs[1])
	}
	if(is.null(df_b)){
		df_b=10
	}
	if(is.null(grid)){
		grid=1000
	}
	##################
	###parameter setting
	#rho=1	### parameter for error process epsilon_t, O-U process
	#tmax= 1001;	### maximum length of time to generated data
	#t=1000	### length of time
	#iter= 100;	### number of replications in the simulation
	#grid=1000;	### the number of grid points used to construct the functional time series X_t and epsilon_t\
	#coef_b=6	### k=6
	#min_obs=40	### the minimal number of observations at each time
	#######################################################K
	x_grid=seq(0,1,by=1/grid)
	coef_b=df_b+1
	b_grid=seq(-1,1,by=1/grid)
	b_grid_sp = ns(b_grid,df=df_b);
	b_grid_full= cbind(rep(1,2*grid+1), b_grid_sp);
	indep= matrix(0, (n+1)*(t-1),coef_b)
	bstar_grid= matrix(0, grid+1, coef_b)
	b_matrix= matrix(0, (grid+1), (grid+1))
	bstar_grid_full= matrix(0, (grid+1)*t,coef_b)
	index=matrix(0,t,n)
	for(i in 1:t){
		for(j in 1:num_obs[i]){
			index[i,j]=which(abs(x_pos[i,j]-x_grid)==min(abs(x_pos[i,j]-x_grid)))
		}
	}
	for (i in 1:(t-1)){
		rec_x= approx(x_pos[i,1:num_obs[i]],f[i,1:num_obs[i]],xout=x_grid,rule=2,method='linear')	
                           ### rec_x is the interpolation of x_t
		for(j in 1:coef_b){
			for (k in 1:(grid+1)){
				b_matrix[k,]= b_grid_full[seq((k+grid),k, by=-1),j]/(grid+1);
			}
			bstar_grid[, j]= b_matrix %*% matrix(rec_x$y,ncol=1,nrow=grid+1)
		}
		bstar_grid_full[(i-1)*(grid+1)+1:(grid+1),]=bstar_grid	###convoluttion of basis spline function and x_t
		tmp=bstar_grid[index[i+1,1:num_obs[i+1]],]
		indep[(i-1)*(n+1)+1:num_obs[i+1],]=tmp
	}


	if(p==1){
		#################AR(1)
		pdf4=function(para4)	###Q function in formula (12)
		{	
			mat1= matrix(0, df_b+1, df_b+1);
			mat2= matrix(0, df_b+1, 1);
			### correlation matrix of error process at observation points
			psi_invall=matrix(0,(n+1)*(t-1),(n+1))
			for(i in 2:t){
				psi= para4;
				tmp.n=num_obs[i]
				psi_mat= matrix(1, tmp.n, tmp.n)
				for(k in 1:(tmp.n)){
					for(j in 1:(tmp.n)){
						psi_mat[k,j]= psi^(abs(x_pos[i,j]-x_pos[i,k]))
					}
				}
				psi_matinv2=solve(diag(weight(x_pos[i,1:tmp.n]))%*%psi_mat%*%diag(weight(x_pos[i,1:tmp.n])))
				psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]=psi_matinv2
				M_t=indep[(i-2)*(n+1)+(1:tmp.n),]
				tmp_mat= t(M_t) %*% psi_matinv2;
				mat1= mat1+ tmp_mat %*% M_t;			### matrix under the first brackets in formula (15)
				mat2= mat2+ tmp_mat %*% f[i,1:tmp.n];	### matrix under the second brackets in formula (15)
			}
			phihat= solve(mat1)%*%mat2
			ehat=0	###e(beta,phi) in formula (12)
			log.mat=0
			for (i in 2:t){
				tmp.n=num_obs[i]
				M_t= indep[(i-2)*(n+1)+(1:tmp.n),]
				eps= f[i,1:tmp.n]- M_t%*% phihat
				psi_matinv= psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]
				epart= t(eps)%*%psi_matinv %*%eps
				ehat= epart+ehat
				log.mat=log.mat+ log(det(psi_matinv))
			}
			l=-sum(num_obs[2:t])/2*log(ehat)+1/2*log.mat	### the number defined in formula (14)
			return(-l);
		}	
		para4= exp(-1)
		result4=optim(para4,pdf4, lower=0.001, upper=0.999, method='L-BFGS-B')
		mat1= matrix(0, df_b+1, df_b+1);
		mat2= matrix(0, df_b+1, 1);
		psi_invall=matrix(0,(n+1)*(t-1),(n+1))
		psi=result4$par
		for(i in 2:t){
			tmp.n=num_obs[i]
			psi_mat= matrix(1, tmp.n, tmp.n)
			for(k in 1:tmp.n){
				for(j in 1:tmp.n){
					psi_mat[k,j]=psi^(abs(x_pos[i,j]-x_pos[i,k]))
				}
			}
			psi_matinv2=solve(diag(weight(x_pos[i,1:tmp.n]))%*%psi_mat%*%diag(weight(x_pos[i,1:tmp.n])))	
			psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]= psi_matinv2
			M_t= indep[(i-2)*(n+1)+(1:tmp.n),]
			tmp_mat= t(M_t) %*% psi_matinv2;
			mat1= mat1+ tmp_mat %*% M_t;
			mat2= mat2+ tmp_mat %*% f[i,1:tmp.n];
		}
		phihat= solve(mat1)%*%mat2	
		phihat_func=matrix(0,p,(1+2*grid))
		for(i in 1:p){
			phihat_func[i,]= cbind(rep(1,2*grid+1),b_grid_sp)%*%phihat[(1+(coef_b*(i-1))):(i*coef_b)];		###b-spline coefficient of estimated phi_1
		}
		ehat=0		###e(beta,phi) in formula (12)
		for (i in 2:t){
			tmp.n=num_obs[i]
			M_t= indep[(i-2)*(n+1)+(1:tmp.n),]
			eps= f[i,1:tmp.n]- M_t%*% phihat
			psi_matinv=psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]
			epart= t(eps)%*%psi_matinv %*%eps
			ehat= epart+ehat
		}
		rho_hat= -log(psi)
		sigma_hat=ehat/sum(num_obs[2:t])*2*rho_hat	

	}
	indep.p=indep	
	if(p>1){
		for(q in 2:p){
			indep.tmp=indep
			for (i in 1:(t-q)){
				for(j in 1:num_obs[i+q]){
					index[i+q,j]= which(abs(x_pos[i+q,j]-x_grid)==min(abs(x_pos[i+q,j]-x_grid)))
				}
				tmp=bstar_grid_full[(i-1)*(grid+1)+index[i+q,1:num_obs[i+q]],]
				indep.tmp[(i-1)*(n+1)+1:num_obs[i+q],]=tmp
			}
			indep.test=cbind(indep.p[(n+2):((n+1)*(t-q+1)),], indep.tmp[1:((n+1)*(t-q)),])
			indep.p=indep.test

		}
	
		pdf4=function(para4)	###Q function in formula (12)
		{	
			mat1= matrix(0, p*(df_b+1), p*(df_b+1));
			mat2= matrix(0, p*(df_b+1), 1);
			### correlation matrix of error process at observation points
			psi_invall=matrix(0,(n+1)*(t-1),(n+1))
			for(i in (p+1):t){
				psi= para4;
				tmp.n=num_obs[i]
				psi_mat= matrix(1, tmp.n, tmp.n)
				for(k in 1:(tmp.n)){
					for(j in 1:(tmp.n)){
						psi_mat[k,j]= psi^(abs(x_pos[i,j]-x_pos[i,k]))
					}
				}
				psi_matinv2=solve(diag(weight(x_pos[i,1:tmp.n]))%*%psi_mat%*%diag(weight(x_pos[i,1:tmp.n])))
				psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]=psi_matinv2
				M_t=indep.p[(i-p-1)*(n+1)+(1:tmp.n),]
				tmp_mat= t(M_t) %*% psi_matinv2;
				mat1= mat1+ tmp_mat %*% M_t;			### matrix under the first brackets in formula (15)
				mat2= mat2+ tmp_mat %*% f[i,1:tmp.n];	### matrix under the second brackets in formula (15)
			}
			phihat= solve(mat1)%*%mat2
			ehat=0	###e(beta,phi) in formula (12)
			log.mat=0
			for (i in (p+1):t){
				tmp.n=num_obs[i]
				M_t= indep.p[(i-p-1)*(n+1)+(1:tmp.n),]	
				eps= f[i,1:tmp.n]- M_t%*% phihat
				psi_matinv= psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]
				epart= t(eps)%*%psi_matinv %*%eps
				ehat= epart+ehat
				log.mat=log.mat+ log(det(psi_matinv))
			}
			l=-sum(num_obs[(p+1):t])/2*log(ehat)+1/2*log.mat	### the number defined in formula (14)
			return(-l);
		}
		para4= exp(-1)
		result4=optim(para4,pdf4, lower=0.001, upper=0.999, method='L-BFGS-B')
		mat1= matrix(0, p*(df_b+1), p*(df_b+1));
		mat2= matrix(0, p*(df_b+1), 1);
		psi_invall=matrix(0,(n+1)*(t-1),(n+1))
		psi=result4$par
		for(i in (p+1):t){
			tmp.n=num_obs[i]
			psi_mat= matrix(1, tmp.n, tmp.n)
			for(k in 1:tmp.n){
				for(j in 1:tmp.n){
					psi_mat[k,j]=psi^(abs(x_pos[i,j]-x_pos[i,k]))
				}
			}
			psi_matinv2=solve(diag(weight(x_pos[i,1:tmp.n]))%*%psi_mat%*%diag(weight(x_pos[i,1:tmp.n])))	
			psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]= psi_matinv2
			M_t= indep.p[(i-p-1)*(n+1)+(1:tmp.n),]
			tmp_mat= t(M_t) %*% psi_matinv2;
			mat1= mat1+ tmp_mat %*% M_t;
			mat2= mat2+ tmp_mat %*% f[i,1:tmp.n];
		}
		phihat= solve(mat1)%*%mat2
		phihat_func=matrix(0,p,(1+2*grid))
		for(i in 1:p){
			phihat_func[i,]= cbind(rep(1,2*grid+1),b_grid_sp)%*%phihat[(1+(coef_b*(i-1))):(i*coef_b)];		###b-spline coefficient of estimated phi_1
		}	
		ehat=0		###e(beta,phi) in formula (12)
		for (i in (p+1):t){
			tmp.n=num_obs[i]
			M_t= indep.p[(i-p-1)*(n+1)+(1:tmp.n),]
			eps= f[i,1:tmp.n]- M_t%*% phihat
			psi_matinv=psi_invall[(i-2)*(n+1)+(1:tmp.n),1:tmp.n]
			epart= t(eps)%*%psi_matinv %*%eps
			ehat= epart+ehat
		}
		rho_hat= -log(psi)
		sigma_hat=ehat/sum(num_obs[(p+1):t])*2*rho_hat
	}

	est_cfarh <- list(phi_coef=t(matrix(phihat,df_b+1,p)),phi_func=phihat_func, rho=rho_hat, sigma2=sigma_hat)
	return(est_cfarh)
}


#' Prediction of CFAR Processes
#'
#' Prediction of CFAR processes.
#' @param model CFAR model.
#' @param f the functional time series data.
#' @param m the forecasting horizon.
#' @references
#' Liu, X., Xiao, H., and Chen, R. (2016) Convolutional autoregressive models for functional time series. \emph{Journal of Econometrics}, 194, 263-282.
#' @return The function returns a prediction of the CFAR process.
#' @examples
#' phi_func= function(x)   
#' {
#'  	return(dnorm(x,mean=0,sd=0.1))
#' }
#' y=g_cfar1(100,5,phi_func)
#' library(MASS)	
#' library(splines)
#' f_grid=y$cfar
#' index=seq(1,1001,by=10)
#' f=f_grid[,index]
#' est=est_cfar(f,1)
#' pred=p_cfar(est,f,1)
#' @export
p_cfar <- function(model, f, m=3){
	###p_cfar: this function gives us the prediction of functional time series. There are 3 input variables:
	###1. model is the estimated model obtained from est_cfar function
	###2. f is the matrix which contains the data
	###3. m is the forecasting horizon

	####################
	###parameter setting
	##t= 100;	### length of time
	##grid=1000;	### the number of grid points used to construct the functional time series X_t and epsilon_t
	#p is the cfar order
	##n=100		### number of observations for each X_t, in the paper it is 'N'
	###############################################################################
	if(!is.matrix(f))f <- as.matrix(f)
	###################
	###Estimation
	n=dim(f)[2]-1
	t=dim(f)[1]
	p=dim(model$phi_coef)[1]
	grid=(dim(model$phi_func)[2]-1)/2
	pred=matrix(0,t+m,grid+1)
	x_grid=seq(0,1,by=1/grid)
	x=seq(0,1,by=1/n)
	for(i in (t-p):t){
		pred[i,]= approx(x,f[i,],xout=x_grid,rule=2,method='linear')$y
	}
	phihat_func=model$phi_func
	for (j in 1:m){
		for (i in 1:(grid+1)){
			pred[j+t,i]= sum(phihat_func[1:p,seq((i+grid),i, by=-1)]*pred[(j+t-1):(j+t-p),])/grid
		}
	}
	pred_cfar=pred[(t+1):(t+m),]
	return(pred_cfar)
}


#' Partial Curve Prediction of CFAR Processes
#'
#' Partial prediction for CFAR processes. t curves are given and we want to predit the curve at time t+1, but we know the first n observations in the curve, to predict the n+1 observation.
#' @param model CFAR model.
#' @param f the functional time series data.
#' @param new.obs the given first n observations.
#' @references
#' Liu, X., Xiao, H., and Chen, R. (2016) Convolutional autoregressive models for functional time series. \emph{Journal of Econometrics}, 194, 263-282.
#' @return The function returns a prediction of the CFAR process.
#' @examples
#' phi_func= function(x)   
#' {
#'  	return(dnorm(x,mean=0,sd=0.1))
#' }
#' y=g_cfar1(100,5,phi_func)
#' library(MASS)	
#' library(splines)
#' f_grid=y$cfar
#' index=seq(1,1001,by=10)
#' f=f_grid[,index]
#' est=est_cfar(f,1)
#' pred=p_cfar_part(est,f,1)
#' @export
p_cfar_part  <- function(model, f, new.obs){
	##p-cfar_part: this function gives us the prediction of x_t(s), s>s_0, give x_1,\ldots, x_{t-1}, and x(s), s <s_0. There are three input variables.
	##1. model is the estimated model obtained from est_cfar function
	##2. f is the matrix which contains the data
	##3. new.obs is x_t(s), s>s_0.
	t=dim(f)[1]
	p=dim(model$phi_coef)[1]
	f_grid=matrix(0,t,grid+1)
	grid=(dim(model$phi_func)[2]-1)/2
	x_grid=seq(0,1,by=1/grid)
	n=dim(f)[2]-1
	index=seq(1,grid+1,by=grid/n)
	x=seq(0,1,by=1/n)
	for(i in (t-p):t){
		f_grid[i,]=approx(x,f[i,],xout=x_grid,rule=2,method='linear')$y
	}
	pred=matrix(0,grid+1,1)
	phihat_func=model$phi_func
	for (i in 1:(grid+1)){
		pred[i]= sum(phihat_func[1:p,seq((i+grid),i, by=-1)]*f_grid[t:(t-p+1),])/grid
	}
	pred_new=matrix(0,n+1,0)
	pred_new[1]=pred[1]
	rho=model$rho
	pred_new[2:(n+1)]=pred[index[1:n]]+(new.obs[1:n]-pred[index[1:n]])*exp(-rho/n)
	return(pred_new)
}
ConvFuncTimeSeries/test3 documentation built on May 29, 2019, 11:41 a.m.