I am trying to simulate a dataset that matches empirical data that I have, but am unsure how to estimate the errors in the original data. The empirical data includes heteroscedacity, but I am not interested in transforming it away, but rather using a linear model with an error term to reproduce simulations of the empirical data.

For example, let's say I have some an empirical dataset and a model:

    n=rep(1:100,2)
    a=0
    b = 1
    sigma2 = n^1.3
    eps = rnorm(n,mean=0,sd=sqrt(sigma2))
    y=a+b*n + eps
    mod <- lm(y ~ n)

using `plot(n,y)` we get the following.
[![enter image description here][1]][1]

However, if I try to simulate the data, `simulate(mod)`, the heteroscedacity is removed and not captured by the model.

I can use a generalized least squares model

    VMat <- varFixed(~n)
    mod2 = gls(y ~ n, weights = VMat)

 that provides a better model fit based on AIC, but I don't know how to simulate data using the output. 

My question is, how do I create a model that will allow me to simulate data to match the original, empirical data (n and y above). Specifically, I need a way to estimate sigma2, the error, using either using a model? 


  [1]: https://i.sstatic.net/QNfBq.png