How to find parameters of AR(1) model with Student-$t$ error terms? I am trying to find parameters of an AR(1) process with error terms $\epsilon\sim$ $t$ distribution. 
With $\{X_1, X_2,...,X_t\}$ given, what I have done so far is to fit the linear regression 
$$
X_t = \beta_0 + \beta_1 X_{t-1},
$$
and then calculate the residuals $\{r_1,r_2,...,r_{t-1}\}$. The model assumption is that $\{r_t\}$ follows a Student-$t$ distribution, however, the standard deviation of $\{r_t\}$  is smaller than 1 (say 0.3). Therefore it could be a scaled $t$ distribution. I was wondering if there is any method for me to find the scaling parameter and the degrees of freedom of the distribution of these residuals?
Also, is it OK to directly use least square method to find $\beta_0$ and $\beta_1$ first? I know they should still be unbiased but no longer the MLE. I was wondering if I could fit the AR(1) model directly with error terms follow a $t$ distribution?
 A: 
Also, is it OK to directly use least square method to find $\beta_0$ and $\beta_1$ first? I know they should still be unbiased but no longer the MLE.

In autoregressive models the OLS parameter estimates are biased (but still consistent).

I was wondering if I could fit the AR(1) model directly with error terms follow a t distribution?

You may use maximum likelihood estimation. Formulate the model likelihood for the scaled Student-$t$ case and find the parameter values that maximize it using optimization. 
(This works fine in theory, but I think I have read somewhere that the estimates of the degrees of freedom can be quite imprecise in practice. I cannot recall the source, though.)
A: You may want to try with vector generalized linear models (VGLMs) applied to time series in R, part of my PhD. My family function AR.studentt.ff(),
through the modelling function vglm(), estimates the three parameters (location, scale and degrees of freedom) of the
the errors distribution, assumed as shift--scaled Student-t, by MLE using Fisher scoring. This family function imposes an AR model on the data, as yours above, also estimated.
My R package, VGAMextra, will be realized soon via CRAN, with further choices. This is an extension of the package VGAM. See below for a very simple example.
Note, only central Student-ts are handled currently.
### Estimate the parameters of the errors distribution for an
## AR(1) model. Sample size = 50

set.seed(20180218)
nn <- 250
y  <- numeric(nn)
ncp   <- 0           # Non--centrality parameter
nu    <- 3.5         # Degrees of freedom.
theta <- 0.45        # AR coefficient
res <- numeric(250)  # Vector of residuals.

y[1] <- rt(1, df = nu, ncp = ncp)
for (ii in 2:nn) {
  res[ii] <- rt(1, df = nu, ncp = ncp)
  y[ii] <- theta * y[ii - 1] + res[ii]
}
# Remove warm up values.
y <- y[-c(1:200)]

### Fitting the model.
AR.stut.er.fit <- vglm(y ~ 1, AR.studentt.ff(order = 1), data = data.frame(y = y), trace = TRUE)

VGLM    linear loop  1 :  loglikelihood = -98.165573
VGLM    linear loop  2 :  loglikelihood = -97.436435
VGLM    linear loop  3 :  loglikelihood = -97.400721
VGLM    linear loop  4 :  loglikelihood = -97.396789
VGLM    linear loop  5 :  loglikelihood = -97.396396
VGLM    linear loop  6 :  loglikelihood = -97.396352
VGLM    linear loop  7 :  loglikelihood = -97.396347
VGLM    linear loop  8 :  loglikelihood = -97.396347

Estimated coefficients of the AR component:
(Intercept)         AR1 
   -0.13496     0.55319 

# Estimated parameters of the Student--t on the errors.
 Coef(AR.stut.er.fit)
 location    scale       df 
 0.053767 1.324088 4.226470 

