Conditionally heteroskedastic linear regression: How can I model variance from given predictors? As a motivating example, consider this hypothesis: A friend and I were discussing how it might be the case that, the more one person is given control in the process of making a movie, the higher the variance of the quality of that film. The more hands that go into making a film, the safer of a bet it is—but the less likely it will be that it is great.
Here are some simulated data that could characterize this hypothesis:
set.seed(1839)
x <- runif(1000, 0, 100)
y <- rnorm(1000, 0, x)
dat <- data.frame(x, y)

library(ggplot2)
ggplot(dat, aes(x = x, y = y)) +
  geom_jitter() +
  theme_light() +
  labs(y = "Quality of Film", 
       x = "Amount of Artistic Control in One Person")


Note that here, I am not interested in predicting the expected value, but the expected variance of $y$ from $x$.
I am familiar with how one can create models for both location and dispersion in beta regression and how one can create models for both location and scale in conditionally heteroskedastic truncated regression.
However, is there a way to model conditional heteroskedasticity in an ordinary least squares regression like this? Are there regression models that exist where one can model the variance of $y$ as a function of $x$? When I try looking for these models, I usually end up with people wanting to use sandwich estimators or weighted least squares or something.
I don't want to treat heteroskedasticity as a problem with the data, but something to be explicitly modeled as a function of predictors $X$. Is there a way to do this with linear regression?
 A: If I understand you correctly, this can be done using the gamlss() function in the gamlss package. 
Given your example, you can model the scale parameter $\sigma$ as follows:
fit <- gamlss(y~x, sigma.formula = ~x, data = dat, family = NO)
> summary(fit)
******************************************************************
Family:  c("NO", "Normal") 

Call:  gamlss(formula = y ~ x, sigma.formula = ~x, family = NO, data = dat) 

Fitting method: RS() 

------------------------------------------------------------------
Mu link function:  identity
Mu Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.32765    1.23130  -0.266    0.790
x           -0.03359    0.04671  -0.719    0.472

------------------------------------------------------------------
Sigma link function:  log
Sigma Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.4663263  0.0533045   46.27   <2e-16 ***
x           0.0254100  0.0009709   26.17   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
No. of observations in the fit:  1000 
Degrees of Freedom for the fit:  4
      Residual Deg. of Freedom:  996 
                      at cycle:  3 

Global Deviance:     10296.29 
            AIC:     10304.29 
            SBC:     10323.92 
******************************************************************

See also Table 1 in the link above for all other available distributions.
A: @Stefan gave rich family of models with the use of gamlss so my response here is more of an extended comment.
For some models with non-constant variance one can just do a simple "rearrangement" of the model and use simpler functions such as lm.
Suppose the model is
$$y=a+b x + x \epsilon$$
with $\epsilon \sim N(0,\sigma^2)$.  We can just divide everything by $x$ and end up with a "standard" linear model with constant variance:
$${y \over x} = {a \over x} + b + \epsilon$$
We now have a typical linear model with $y/x$ as the dependent variable, $1/x$ is the independent variable with "slope" $a$ and intercept $b$ and the variance is constant.
