Is there a "beta distribution" over the entire real line? The malleability of the beta distribution pleases me; it can be symmetric, asymmetric, platykurtic and so on, as the following picture shows us:

I thought it would be interesting to use it to model the errors in regression problems, but its support is only defined to be $(0,1)$ so I would like to know if there is a "beta" distribution with support $(-\infty,\infty)$ because a error lying on $(0,1)$ is such a strong assumption.
 A: It seems to me the convex shape the beta distribution can take is a result of the bounded sample space.  If we move to a sample space that is bounded on one side then we are dealing with distributions like gamma, Weibull, etc.  If we move to an unbounded sample space then we of course have the normal distribution.
All of these distributions, including the beta and scaled beta, can be used in regression via generalized linear models.  Rather than developed around an "error term" these models are typically developed by expressing a parameter in the likelihood or probability distribution in terms of a linear predictor.  You could re-parameterize a model in terms of its mean or median and express this in terms of a linear predictor.
A: One way to make a "real-valued beta distribution" would be to transform the interval $(0,1)$ onto the real line. One way of doing that is the logistic function
$$ \text{logit}(x)= \log(\frac{x}{1-x}) $$ known from, for example, logistic regression. So if $X \sim \mathcal{Beta}(\alpha,\beta)$, let
$$ Y=\text{logit}(X)= \log(\frac{X}{1-X}). $$
Then $Y$  will be a random variable with range on the full real line, with a distribution we can call logistic-Beta.
By standard transformation methods we find the density function of $Y$
as
$$
f_Y(y) =\frac{e^{\alpha y}}{B(\alpha, \beta) (1+e^y)^{\alpha+\beta}}
$$ where $B$ is the beta function. This will give a flexible family of (unimodal) distributions on the real line. In this comprehensive reference that I just found, this distribution is called the Beta-Logistic Distribution. That is a useful search term!
An example is below:

The same transformation could be used with mixtures of betas, for instance, to get even multi-modal distributions. In reality, this example is a special case of the log-ratio transformation used with Compositional data,  see also How to perform isometric log-ratio transformation.
Transforming a beta random variable to the positive line (in a specific way) only gives the beta prime distribution, a recent post with an example is Distribution of the exponential of an exponentially distributed random variable?.
Code use for the figure is below:
dlogisticBeta <- function(x, alpha, beta, log=FALSE) {
    stopifnot( (alpha>0)&&(beta>0) )
    logans <- alpha*x - lbeta(alpha,  beta) -
        (alpha + beta)*log1p( exp(x))
    if(log) return(logans) else return(exp(logans))
}

oldpar <- par(mfrow=c(1, 2))
plot( function(x) dlogisticBeta(x, 0.5, 0.5), from=-10,  to=10, col="red", 
     main= expression(paste(alpha == 0.5,",   ",   beta==0.5)) )
plot( function(x) dlogisticBeta(x, 0.5, 3), from=-10,  to=10, col="red",  
     main= expression(paste(alpha == 0.5,",   ",   beta==3)) )
par(oldpar)

A: Beta regression is used to model continuous proportions; & that such observations are bounded by 0 & 1 is as reasonable an assumption as it is a strong one (if you'd call it an assumption rather than a definition). As noted in @GeoffreyJohnson's answer, it's the conditional distribution of each observation that's taken to be a beta distribution, not the error.  In a typical beta regression you reparametrize the beta distribution with mean $\mu = \frac{\alpha}{\alpha + \beta}$ & precision $\nu = \alpha + \beta$; the joint density $f_\boldsymbol{Y}$ of $n$ independent observations $\boldsymbol{y}$ is given by $$f_\boldsymbol{Y}(\boldsymbol{y})=\prod_i^n\frac{\Gamma(\nu)}{\Gamma(\mu_i\nu)\Gamma\left[(1-\mu_i)\nu\right]} \cdot y_i^{\mu_i\nu-1}(1-y_i)^{(1-\mu_i)\nu-1}$$
The link function $g$ connects the predictor values $\boldsymbol{x}_i$ & the coefficients $\boldsymbol\beta_i$ to the mean for each observation $\mu_i$:
$$g(\mu_i) = \eta(\boldsymbol{x}_i,\boldsymbol\beta_i)=\boldsymbol{x}_i^\mathrm{T}\boldsymbol{\beta}$$
The precision is common across all observations.
Each row of plots below shows therefore how the distribution may be changed by changing the predictor values:

Now there's evident flexibility ('malleability') there, but there are also evident constraints, among them these:

*

*$f_Y(y; \mu,\nu)=f_Y(1-y; 1-\mu,\nu)$. This provides a desirable equivariance—if Dr A models the proportion of the day an animal spends asleep & Dr B the proportion it spends awake, we'd hope they'd make substantively identical inferences/predictions.


*$\operatorname{Var} Y_i =\frac{\mu_i(1-\mu_i)}{1+\nu}$. The variance is at a maximum when $\mu=\frac{1}{2}$, increasing as $\mu\rightarrow 0$ or $\mu\rightarrow 1$. Qualitatively this is what you'd often expect from a continuous proportion, squashing up against the bounds.


*$\frac{\operatorname{E} (Y_i - \operatorname{E} Y_i)^3}{\left(\sqrt{\operatorname{Var} Y}\right)^3}  =\frac{2(1-2\mu_i)\sqrt{1+\nu}}{(2+\nu)\sqrt{\mu(1-\mu)}}$. The skewness is positive when $\mu<\frac{1}{2}$, increasing as $\mu\rightarrow$; it's negative when $\mu>\frac{1}{2}$, decreasing as $\mu\rightarrow 1$; & of course it's zero when $\mu=\frac{1}{2}$. Again, this is qualitatively what you'd often expect.
The point of going into all this is first that in a sense there's no free lunch—you pay (in free parameters) for a certain amount of flexibility, & allocating it to some distributional properties means constraining others—and second that it's not at all clear what you'd want to carry over to modelling a real-valued response, beyond allowing skew & kurtosis to vary in some fashion (as @Glen_b comments).
