# Function to compute the standard errors assuming a
# Gaussian noise process
## NOT WORKING YET

armavar<-function(n,kar,kma,sigma2){
	p<-length(kar)
	q<-length(kma)
	k<-p+q
	Gammapq<-matrix(0,k,k)
	sigma<-sqrt(sigma2)
# Include different noise variance in AR, MA model specification
	x<-arima.sim(model=list(ar=kar),10000,rnorm(n=10000,sd=sigma))
	y<-arima.sim(model=list(ma=kma),10000,rnorm(n=10000,sd=sigma))
	aracf<-acf(x,plot=F,lag=p,type="covariance")
	maacf<-acf(y,plot=F,lag=q,type="covariance")
	armaacf<-acf(cbind(x,y),plot=F,lag=max(q,p),type="covariance")
# Fill in AR part	
	for (i in 1:p){
		for (j in 1:p) Gammapq[i,j]<-aracf$acf[abs(i-j)]
		}
# Fill in MA part
	for (i in (p+1):k){
		for (j in (p+1):k) Gammapq[i,j]<-maacf$acf[abs(i-j)]
	}
# Fill in ARMA part
	for (i in 1:p){
		for (j in (p+1):k) Gammapq[i,j]<-armaacf$acf[abs(i-j),1,1]
	}	
	Gammapq[(p+1):k,1:p]<-t(Gammapq[1:p,(p+1):k])
Betavar<-sigma2*solve(Gammapq)/n
return(Betavar)
}