Why gls() is giving the same estimates for the same dataset independently of the input arguments?

I'm currently struggling with linear autocorrelated models. I have correctly simulated some datasets to understand how homoscedasticity and heteroscedasticity work and now I'm focused on correlated residuals. Therefore I do simulate a dataset with correlated residuals.

The problem: in the end I compare the estimates of the models that do not or do include a specification to model the correlation structure, expecting the latter one to obtain better results and/or estimates. However, all the estimates are the same (?)

My simulated dataset:

• 5 subjects
• 3 days
• correlation structure = compound symmetry. Rho value = 0.85 (I've also tried AR1, and happens the same).
• Variance function: the variance is depending on the mean. Therefore the model is heterocedastic. More specifically, the variance of the residuals = mean of the y values.

Code:

# DATA SIMULATION
days <- gl(3,1,15) # 3 days.
subject <- gl(5,3,15) # subjects
dm <- model.matrix(~ 1+ days, data.frame(days))  # design matrix
beta <- c(10,20,30)
mu <- dm%*%beta
# variance-covariance matrix
CS <- corCompSymm(value = 0.85, form = ~ days|subject, fixed = FALSE)
CS <- Initialize(CS, data =data.frame(days))
corMatrixCS<- bdiag(corMatrix(CS))
sd <- diag(sqrt(as.vector(mu)))

# residuals
set.seed(1234)
(eCSV <- mvrnorm(1, mu = rep(0,15), Sigma = sd %*% as.matrix(corMatrixCS) %*% sd))
yCSV <- ( dm %*% beta ) + eCSV


Plot of the simulated data (y values vs days).

Input arguments of the gls function.

mB <- gls( model = y ~ 1, data = dataCorrVar) # offset model
mB2 <- gls( model = y ~ 1 + days, data = dataCorrVar) # No correlation or variance applied
mB3 <- gls(model =  y ~ 1 + days, data = dataCorrVar,
correlation = corCompSymm(value = 0.85, form=~ days|subject) )
mB4 <- gls( model = y ~ 1 + days, data = dataCorrVar, weights = varPower(fixed = 1))
mB5 <- gls(model =  y ~ 1 + days, data = dataCorrVar,
correlation = corCompSymm(value = 0.85, form=~ days|subject),
weights = varPower(fixed = 1))


Output table:

1st model: offset model. Only the intercept is calculated.

2nd model: "raw model". No variance or correlation modeling.

3rd model: Only correlation is modeled.

4th model: Only variance is modeled.

5th model: Both variance and correlation structure are modeled.

As I said before, I expect the 4th model to give the best results. And it does in terms of AIC or BIC values. However, to me, it is extremely rare that all the estimates are exactly the same and only the standard errors are modified (std error between parenthesis).