Do parameters stay unchanged when GARCH is scaled? Let's say we have a GARCH($1,1$) process specified as follows:
$y_t = \epsilon_t \sqrt h_t, \quad \epsilon_t \sim N(0,1) \quad \text{i.i.d.}$
$h_t = a_0 + a_1 y^2_{t-1} + b_1 h_{t-1}.$
If we were trying to estimate the parameters $\Theta = (a_0, a_1, b_1)$, and we have a sample $\{y_k\}^n_{k=1}$, would it make a difference if we estimated the parameters from $\{y_k\}^n_{k=1}$ or $100\times\{y_k\}^n_{k=1}$? I basically want to know which (if any) parameters would also scale. I am rescaling the $y$ for numerical stability.    
 A: $a_0$ will scale quadratically with the scaling of $y_t$, while $a_1$ and $b_1$ will stay the same. Here is why. Take the original model
\begin{aligned}
y_t &= \sqrt h_t \epsilon_t, \\
h_t &= a_0 + a_1 y^2_{t-1} + b_1 h_{t-1}, \\
\epsilon_t &\sim i.i.N(0,1) \\
\end{aligned}
and scale $y_t$ with a positive constant $c$. This turns the model into
\begin{aligned}
cy_t &= \sqrt{c^2 h_t} \epsilon_t, \\
c^2 h_t &= c^2 a_0 + a_1 (cy_{t-1})^2 + b_1 (c^2 h_{t-1}), \\
\epsilon_t &\sim i.i.N(0,1) \\
\end{aligned}
which can be expressed as
\begin{aligned}
\tilde y_t &= \sqrt{\tilde  h_t} \epsilon_t, \\
\tilde  h_t &= \tilde  a_0 + a_1\tilde y^2_{t-1} + b_1 (\tilde  h_{t-1}), \\
\epsilon_t &\sim i.i.N(0,1) \\
\end{aligned}
where $\tilde y_t = c y_t$ and $\tilde a_0 = c^2 a_0$.

However, numeric stability issues may create an impression that this is not the case. Try the following code (illustrated with two different R packages for GARCH models) and note how the estimated parameters do not scale all that nicely. (This might be due to the fact that the randomly generated data is i.i.d., not GARCH, making estimated the GARCH coefficients unstable.)
n=1000
set.seed(1); x=rnorm(n) # generate i.i.d N(0,1) data

# With ARMA(1,1) conditional mean
library(rugarch)
spec=ugarchspec()
fit2 =ugarchfit(data=x*100,spec=spec); coef(fit2 ) # scaled by 100
fit1 =ugarchfit(data=x*10 ,spec=spec); coef(fit1 ) # scaled by 10
fit0 =ugarchfit(data=x    ,spec=spec); coef(fit0 ) # raw data
fit_1=ugarchfit(data=x/10 ,spec=spec); coef(fit_1) # scaled by 0.1
fit_2=ugarchfit(data=x/100,spec=spec); coef(fit_2) # scaled by 0.01

# With constant conditional mean
library(garchx)
fit2 =garchx(y=x*100); coef(fit2 )
fit1 =garchx(y=x*10 ); coef(fit1 )
fit0 =garchx(y=x    ); coef(fit0 )
fit_1=garchx(y=x/10 ); coef(fit_1)
fit_2=garchx(y=x/100); coef(fit_2)
```

A: Define $\tilde{y}_t = cy_t$ with $c > 0$ and you get 
$$
\tilde{y}_t = \sqrt{\tilde{h}_t}\epsilon_t
$$
where $\tilde{h}_t = c^2 h_t$ is the new volatility process. The dynamics for this scaled up volatility process can be obtained by multiplying both sides of the original equation by $c^2$:
\begin{align*}
\tilde{h}_t &= c^2 a_0 + c^2 a_1 y^2_{t-1} + c^2 b_1 h_{t-1} \\
&= \tilde{a}_0 + a_1 \tilde{y}^2_{t-1} + b_1 \tilde{h}_{t-1}.
\end{align*}
The volatility is a linear function of its past value, and the past squared observation. You can see that $a_1$ and $b_1$ stay the same, but $a_0$ turns into $\tilde{a}_0 = c^2 a_0$. 
