Using a 2nd neural network to predict 1st neural network prediction error So for example, we are trying to predict the amount of rainfall in the afternoon base on continuous features such as humidity and temperature in the morning.

1st neural network: Regression neural network on features to give one output label which is a continuous value for predicted rainfall in the afternoon.

2nd neural network: Features will be the same as the 1st neural network. But now, the label will be the absolute difference between the 1st neural network predicted rainfall and the actual rainfall.

From this, we can train the 2nd neural network to recognise what particular set of features will result in the 1st neural network giving a 'bad' prediction and be more wary of that 'bad' prediction. In a way, this is like using another neural network to give the confidence level of the first neural network, solely based on the same features (the humidity and temperature in the morning).

I could not find much literature on this subject and am wondering if this idea makes sense in the first place? Perhaps stacking neural networks over each other is a bad idea because it compounds the error from one network to another?

I tried this with some data except my 2nd neural network is a classifier which classifies if the error is above a certain threshold (bad prediction) or below a certain threshold (good prediction).

However, from a few different model runs, it seems that my 2nd neural network usually gives a matthew's correlation coefficient of about 0. This means my 2nd neural network is as good as guessing whether the 1st neural network prediction is good or bad.

So I am not sure if the problem is the idea itself or that my model hyperparameters are bad. 

More details: I used 10 fold cross-validation for the 1st model to get a predicted rainfall for all the data. Then I used another separate 10 fold cross-validation and a siamese neural network for the 2nd model to predict whether the 1st neural network prediction is good or bad.
It's not an unreasonable approach, but I suspect you'd need to define a custom loss function to make it work well.

I can also suggest two different, more sophisticated approaches: the bootstrap, or a variational neural network.

In the bootstrap, you train many classifiers (say, 100 of them); each is trained on a different random sample of the training set, and then you look at the distribution of outputs from these classifiers when you feed in the input $x$ to each of them.

In a variational network, instead of outputting a single number $y$ for the prediction in response to the input $x$, the network outputs two parameters $(\mu,\sigma)$, with the idea that the network is predicting a Gaussian distribution $\mathcal{N}(\mu,\sigma^2)$ as an approximation for $p(y|x)$.  Then you can use this to get a sort of confidence interval for the prediction, e.g., $[\mu-2\sigma,\mu+\sigma]$.  I'm not an expert on this, but I think a variational network is actually very close to what you suggested; we can think of it as two networks, one that outputs $\mu$ (your first network) and one that outputs $\sigma$ (your second network).  However, variational networks are trained with a special loss function.  In particular, if we have an instance $(x_i,y_i)$ in the training set, the loss for the network is the log likelihood $-\log p(y_i)$ where here $p(y_i)$ represents the probability of getting the output $y_i$ from a Gaussian distribution with parameters $\mathcal{N}(\mu,\sigma^2)$, where $\mu,\sigma$ are the two outputs from the network.  Since the Gaussian distribution has probability density function

$$p(y_i) = c e^{-(y_i-\mu)/2\sigma}$$

where $c$ is a constant.  Therefore, the loss function is

$$-\log p(y_i) = (y_i-\mu)/2\sigma + c'$$

where $c'$ is a constant that can be ignored.  Thus, I think you can think of the variational approach as being equivalent to your approach, but with a custom loss function chosen to be appropriate for what you're trying to achieve.