inno.ma<-function(q,estimated.gamma,nmax)
{
  theta=array(NA,nmax)
  t=array(NA,c(nmax,q))
  v=array(NA,nmax+1)

  v[1]=estimated.gamma[1]
  for(n in 1:nmax)
    {
      for(k in max(n-q,0):(n-1))
		 {
		   t[n,n-k]=estimated.gamma[n-k+1]
		   if(k>0 && k-1>n-q)
		     {
		       for(j in max(n-q,0):(k-1))
		 		 t[n,n-k]=t[n,n-k]-t[k,k-j]*t[n,n-j]*v[j+1]
		     }
		   t[n,n-k]=t[n,n-k]/v[k+1]
		 }
      
      v[n+1]=estimated.gamma[1]
      
      for(j in max(n-q,0):(n-1))
		 v[n+1]=v[n+1]-t[n,n-j]^2*v[j+1]
      
      if(v[n+1]<0.0)
		 {
		   for(j in 1:q)
		     theta[j]=NA
		 }
    }

  for(j in 1:q)
    theta[j]=t[nmax,j]
  sig2=v[nmax+1]

  theta
}