I am new to Stan, and I am trying to fit a CFA model with multivariate normal indicator variables. I am using Stan directly instead of blavaan because I am going to build on this model to include a hierarchical part to it. So the CFA part is just a start.

So these are the dimensions of my data.
Please share your Stan program and accompanying data if possible.

When including Stan code in your post it really helps if you make it as readable as possible by using Stan code chunks (```stan) with clear spacing and indentation. For example, use

data {
int N; //number of persons
int ind; //number of indicator variables i.e. 9 with 3 variables per factor
matrix[N, ind] Y; // my input data
int Nf; //number of factors which will be 3
}
parameters {
matrix[2, Nf] ld; //only 2 loadings per factor since a reference variable will be defined for each factor
vector[ind] intercepts;
cholesky_factor_corr[Nf] Nf_corr;
vector < lower = 0 > [ind] sd_res; //residual sds
}
transformed parameters {
matrix[N, ind] mu_Y; //expected values of Y
matrix[N, Nf] LVs; //latent variables
vector[Nf] zeros;
cov_matrix[ind] Sigma;
zeros = rep_vector(0, Nf)
}
model {
//priors
for (k in 1:ind) {
intercepts[k] ~ normal(0, 10);
sd_res [k] ~ cauchy(0, 2.5);
}
for ( i in 1:2) {
for(j in 1:Nf) {
ld[i,j] ~ normal(0,10);
}
}
// distribution of data
to_vector(Y) ~ multi_normal(to_vector(mu_Y), Sigma);
// distribution of factors
to_vector(LVs) ~ multi_normal(zeros, Nf_corr);
}

I keep getting errors based on the multi_normal_lpdf such as: Size of random variable (2709) and rows of covariance parameter (9) must match in size.
I have also left out some part of the code just to keep it simple.
I will be appreciate any help offered or if any one can point me in the direction of a STAN model in hierarchical SEM.

It may help you to look at page 34 of the Stan User’s Guide for version 2.21. There is a model there that uses a multivariate normal, though for continuous predictors.

data {
int<lower=1> K;
int<lower=1> J;
int<lower=0> N;
vector[J] x[N];
vector[K] y[N];
}
parameters {
matrix[K, J] beta;
cov_matrix[K] Sigma;
}
model {
vector[K] mu[N];
for (n in 1:N)
mu[n] = beta*x[n];
y ~ multi_normal(mu, Sigma);
}

In that case there are N observations of y and x. Each y is a vector of K elements and so y is declared as

vector[K] y[N];

Each x is a vector of J predictors and x is declared as

vector[J] x[N];

The coefficients beta that are applied to x to predict the mean y values is a K x J matrix

matrix[K,J] beta;

The predicted mean values, mu, have the same kind of declaration as y

vector[K] mu[N];

A for loop is used to generate the N vectors of mu and then the y values are predicted with

y ~ multi_normal(mu, Sigma);

where Sigma is a K x K covariance matrix. That is vectorized and does not show explicitly that each y[n] which is a K vector is being predicted using a K vector of mu values. In your code, you have changed the Y and mu values into single vectors, which is not what multi-normal is expecting. Because your Y has 2709 values, it expects Sigma to be a 2709 x 2709 matrix. You need your Y to be an array of N vectors each with 9 values.