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.

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).
