Opposite results using Bayesian (STAN) vs Multilevel model (nlme). How is this possible?

Question

My datasets contains the median wages and the cumulative installed wind-capacity for 4000 counties over a period of 20 years. The wages tend to rise over the period and the capacity tends to highly differ between counties. As a first step I want to analyse the broad effect and later group counties by states and other categories.

I ran base level analysis in R using Bayesian Markov Chain Monte-Carlo (MCMC) simulation with STAN and got a small positive correlation which is also backed by other studies. The code must be correct as it is a near 100% replica of another study.

Now i tried to replicate the same using Hierarchical Linear Models using nlme and got a totally opposite result.

nmle:

 model <- lme(median_wage ~ cum_installed_wind + period,
                    random = ~ period | county_id,
                    correlation = corAR1(form = ~ period | county_id),
                    data = obs_level,
                    method = "ML")

Stan:

data {
  int<lower = 0> N; // number of observations
  int<lower = 0> C; // number of counties
  int<lower = 1, upper=C> county[N]; // county ID
  vector[N] median_wage; // median yearly wages
  vector[N] period; // time period
  vector[N] capacity; // cumulative wind capacity
}

transformed data {
  vector[N] st_median_wage;
  vector[N] st_capacity;
  vector[N] st_period;

  // Standardize inputs
  st_median_wage = (median_wage - mean(median_wage)) / sd(median_wage);
  st_capacity = (capacity - mean(capacity)) / sd(capacity);
  st_period = (period - mean(period)) / sd(period);
}

parameters {
  real<lower=0> sigma_alphaW; // random intercept std dev
  real<lower=0> sigma_beta0W; // random slope std dev
  real<lower=0> sigma_y; // residual std dev

  real mu_alphaW; // mean of intercept
  real mu_beta0W; // mean of slope for period
  real beta1W; // fixed effect of wind capacity

  vector[C] alphaW_raw; // raw random intercepts
  vector[C] beta0W_raw; // raw random slopes for period
}

transformed parameters {
  vector[C] alphaW;
  vector[C] beta0W;
  vector[N] y_hat;

  // Compute random effects
  alphaW = mu_alphaW + sigma_alphaW * alphaW_raw;
  beta0W = mu_beta0W + sigma_beta0W * beta0W_raw;

  // Define the model prediction
  for (i in 1:N) {
    y_hat[i] = alphaW[county[i]] + beta0W[county[i]] * st_period[i] + beta1W * st_capacity[i];
  }
}

model {
  // Priors
  mu_alphaW ~ normal(0, 5);
  mu_beta0W ~ normal(0, 5);
  beta1W ~ normal(0, 5);
  sigma_alphaW ~ cauchy(0, 5);
  sigma_beta0W ~ cauchy(0, 5);
  sigma_y ~ cauchy(0, 5);

  alphaW_raw ~ normal(0, 1);
  beta0W_raw ~ normal(0, 1);

  // Likelihood
  median_wage ~ normal(y_hat, sigma_y);
}

1. Is this code correct for my analysis?
2. How are these opposite results possible?