The main thrust of the article is the development of a joint stochastic model for electricity demand, and wind and solar power production in a given region. The model hinges on special statistical data analysis techniques including the estimation of heavy tail distributions, graphical LASSO fitting procedures, and conditional Monte Carlo simulations. Assuming the availability of point forecasts, we model the deviations from these forecasts instead of modeling the actual quantities of interest. The resolution of the model is determined by the resolution of the forecast data. For the sake of illustration, we implement our model and the corresponding simulation algorithms on data made available by NREL for the Texas region with hourly time resolution, load data at the zone level, and wind and solar power production at the generation asset level. Our numerical simulations confirm that the dependencies identified through the fitting algorithm are consistent with the relative locations of the production assets and the geographical load zones over which the data were collected.