I think we had it right the first time. error_s is a just a vector of i.i.d standard normals
error_s = pm.Normal('error_s',mu=0.0,sd=1, shape=len(returns))
Though I’m still getting an error about input dimension
ValueError: Input dimension mis-match. (input[0].shape[0] = 3812, input[1].shape[0] = 3813)