@model function SSM(n, x0, A, B, Q, P) 
 
 	 # x is a sequence of hidden states 
 	 x = randomvar(n) 
 
 	 # y is a sequence of clamped observations 
 	 y = datavar(Vector{Float64}, n) 
 
 	 # `~` expression creates a probabilistic relationship
 	 # between random variables
 	 x_prior ~ MvNormal(μ = mean(x0), Σ = cov(x0)) 
 	 x_prev = x_prior 
 
 	 # Build the state-space model 
 	 for i in 1:n 
 		   x[i] ~ MvNormal(μ = A * x_prev, Σ = Q) 
 		   y[i] ~ MvNormal(μ = B * x[i], Σ = P) 
 		   x_prev = x[i] 
 	 end 
 end 
 
 observations = load_dataset() 
 
 result = inference(
 	model = SSM(length(observations), x0, A, B, Q, P), 
 	data  = (y = observations,)
 )