```
data {
int<lower=0> N;
int<lower=0> K;
matrix[N,K] y;
matrix[N,K] X;
}
parameters {
cov_matrix[K] lambda_hat;
cov_matrix[K] sigma_0;
cov_matrix[K] L;
vector[K] a;
}
transformed parameters {
matrix[K,N] mu;
vector[4] v;
for (i in 1:K)
v[i]=1;
for (i in 1:N)
mu[,i]=lambda_hat*X[i,]';
}
model {
//sigma_0~ inv_wishart(K,diag_matrix(v));
lambda_hat~ wishart(K,diag_matrix(v));
a~ multi_normal(v,diag_matrix(v));
for (t in 100:N)
y[t,]'~ multi_normal(a+mu[,t],diag_matrix(v));
}
```

In each iteration in t, i.e.,

```
for (t in 100:N)
y[t,]'~ multi_normal(a+mu[,t],diag_matrix(v));
```

I want to compute the variance covariance matrix of y[t,]β, Var(y[1:t,]β)=H(t), where H(t) is data-dependent (data on y upto time t), is there any way of computing each H(t) based on data.

Please help me out.

My R code corresponding to data is,

```
data1<-data[order(nrow(data):1),]
prices<-data1[,c(5,7,9,11)]
orderflow<-data1[,c(6,8,10,12)]
dataList = list(X=as.matrix(orderflow), y=as.matrix(prices), N=48614, K=ncol(orderflow))
```