# Note that n is the sample size and n0 is the burnin

# Uses normal innovations

bilinear <- function(phi,b,n,n0=500) {
	y <- rep((n+n0))
	w <- rnorm(n + n0)
        
	for (t in 2:(n+n0)){
            
		y[t] <- phi * y[t-1] + b * w[t-1] * y[t-1] + w[t]

            }
        
         y = y[-c(1:n0)]

        return(y)
}