Is it possible to generate numbers from Gamma distribution (with parameters shape=10, scale=15, say) which also follow a AR(1) process, simultaneously? If it's possible, than how to do that?
I don't think this sort of thing will generally be possible without changing some aspects of the usual definition of an AR.
[most of original answer removed, since I have better information]
Richard Hardy points out in the comments below that probably meant a Gamma marginal distribution (... and rereading your question, that makes sense); Richard thought to use an innovation-approach to create the model. While it's not possible to pick just any distribution for innovations, there is at least one such innvation-model that does work (see the Lawrance reference below).
Here's some references that will prove useful:
Grunwald G.K., Hyndman R.J., Tedesco, L. and Tweedie, R. L. (2000)
"Non-Gaussian Conditional Linear AR(1) Models,"
Australian & New Zealand Journal of Statistics, 42:4 (Dec) p479-495 DOI: 10.1111/1467-842X.00143 (Working paper here: http://robjhyndman.com/papers/clar.pdf)
This and the research report below discuss a variety of approaches to non-Gaussian AR(1) models. [The above paper gives a general formulation of non-Gaussian AR(1) models that includes nearly all published non-Gaussian AR(1) models, while the report below is partly a survey paper as well as beginning the synthesis of the paper above.]
Grunwald, G.K., Hyndman, R.J. and Tedesco, L.M. (1995),
"A unified view of linear AR(1) models,"
Research Report, Department of Statistics, University of Melbourne
Lawrance, A.J. (1982),
"The innovation distribution of a gamma distributed autoregressive process,"
Scand. J. Statist., 9, 234–236
[You could perhaps more readily construct an AR in the log of a gamma. Fitted to log-data, that would correspond to a shifting scale in the conditional distribution for the gamma.]