Can PyStan models be used to make predictions after fitting? Once I've fit a PyStan model, how can I use it to make future predictions without re-fitting?
 A: This question is really about Stan rather than pystan, so I will demonstrate how one might be able to use the posterior samples to make predictions.
Say you have a very simple Bayesian model for linear regression.
$$ y \sim \mathcal{N} (\alpha +X \beta, \sigma^2) $$
$$ \beta \sim \mathcal{N}(0, 1) $$
$$ \sigma \sim \operatorname{Half-Cauchy}(0, 1) $$
A Stan model for this data might be...
data{
    int n;
    int p;
    matrix[n,p] X;
    vector[n] y;
}
parameters{
    vector[p] beta;
    real alpha;
    real<lower=0> sigma;
}
model{
    beta ~ std_normal();
    sigma ~ cauchy(0,1);
    y ~ normal_id_glm(X, alpha, beta, sigma);
}

We can fit this model on some data.  Again, keep in mind I'm using cmdstanpy not pystan so the setup and execution are not exactly the same, but the ideas still apply.

with open('tmp_model.stan', 'w') as file:
    file.write(stan_model_code)
    
stan_model = cmdstanpy.CmdStanModel(stan_file='tmp_model.stan')

n, p = 100, 5
X = np.random.normal(size = (n,p))
beta = np.ones(p)
y = X@beta + np.random.normal(size = n)

model_data = dict(n=n, p=p, X=X, y=y)

fit = stan_model.sample(model_data)

The fit object contains posterior samples.  If I want to make predictions for a new observation, I need to apply the samples to the new observation, like so
#New observation
new_x = np.zeros(p)

# Extract samples as a dictionary
samples = fit.stan_variables()

# Extract samples for the intercept and coefficents of the linear model
intercept_samples = samples['alpha']
coef_samples = samples['beta']

# Initialize a place to store the predictions
predictions = np.zeros_like(intercept_samples)

# Loop over each sample
for i, (a, b) in enumerate(zip(intercept_samples, coef_samples)):
    
    # Compute the prediction manually
    predictions[i] = a + new_x@b

You can now take the mean, quantiles, etc on the prediction object.

If you're using cmdstanpy, this is easier with the generated quantities block

stan_model_code = '''
data{
    int n;
    int p;
    matrix[n,p] X;
    vector[n] y;
    
    int npred;
    matrix[npred, p] Xpred;
}
parameters{
    vector[p] beta;
    real alpha;
    real<lower=0> sigma;
}
model{
    beta ~ std_normal();
    sigma ~ cauchy(0,1);
    y ~ normal_id_glm(X, alpha, beta, sigma);
}
generated quantities{

    vector[npred] ypred = alpha + Xpred*beta;

}

'''

with open('tmp_model.stan', 'w') as file:
    file.write(stan_model_code)
    
stan_model = cmdstanpy.CmdStanModel(stan_file='tmp_model.stan')

n, p = 100, 5
X = np.random.normal(size = (n,p))
beta = np.ones(p)
y = X@beta + np.random.normal(size = n)
# Place holders
npred = 1
Xpred = np.zeros((1,p))
model_data = dict(n=n, p=p, X=X, y=y, npred=npred, Xpred=Xpred)

fit = stan_model.sample(model_data)


pred_data = model_data.copy()
pred_data['npred'] = 10 # 10 new observations to predict
pred_data['Xpred'] = np.random.normal(size = (10, p)) # New observations
predictions = stan_model.generate_quantities(pred_data, fit)

I'm sure pystan has something similar, but I stopped using pystan a while ago in favor of cmstanpy.
