I am trying to build a simple VAR(p) model using pymc3, but I'm getting some cryptic errors about incompatible dimensions. I suspect the issue is that I'm not properly generating random matrices. Here is an attempt at VAR(1), any help would be welcome:
# generate some data
y_full = numpy.zeros((2,100))
t = numpy.linspace(0,2*numpy.pi,100)
y_full[0,:] = numpy.cos(5*t)+numpy.random.randn(100)*0.02
y_full[1,:] = numpy.sin(6*t)+numpy.random.randn(100)*0.01
y_obs = y_full[:,1:]
y_lag = y_full[:,:-1]
with pymc3.Model() as model:
beta= pymc3.MvNormal('beta',mu=numpy.ones((4)),cov=numpy.ones((4,4)),shape=(4))
mu = pymc3.Deterministic('mu',beta.reshape((2,2)).dot(y_lag))
y = pymc3.MvNormal('y',mu=mu,cov=numpy.eye(2),observed=y_obs)
The last line should be
y = pm.MvNormal('y',mu=mu.T, cov=np.eye(2),observed=y_obs.T)
MvNormal interprets the last dimension as the mvnormal vectors. This is because the behaviour of numpy indexing implies that y_obs is a vector of length 2 containing vectors of length 100 (y_lag[i].shape == (100,))