Here is an argument in two steps. Firstly, $M = pR + (1-p)Q$ is the variance-covariance matrix of some vector of random variables. Proof: we require to show that $M$ is semi-positive definite. We know that $R$ and $Q$ are, since they are correlation matrices. Hence for any non-zero vector $x$, $x^tPx \geq 0$ and $x^tQx \geq 0$. Since $p\geq0$ and $(1-p)\geq0$ we have, as required:
$$x^tMx = p(x^tRx) + (1-p)(x^tQx) \geq 0$$.
For this covariance matrix to be a correlation matrix we additionally require that the random variables that form the components of the vector have variance one. The variances are the diagonal elements, and we already know that the diagonal elements of $R$ and $Q$ are one since they are correlation matrices. So:
$$M_{ii} = pR_{ii} + (1-p)Q_{ii} = p(1) + (1-p)(1) = 1$$
This estalishes that "the weighted average of two correlation matrices is a correlation matrix", which is the question in the title.
It doesn't address the question implicit in the body of the post though:
I claim that $\sqrt pX+\sqrt{1-p}Y$ is a random vector that has
correlation matrix $pR+(1-p)Q$. Here, $X$ is the random vector with
correlation $R$ and $Y$ is the random vector with correlation $Q$.
Also, $X$ and $Y$ are independent. I can get the covariance part to
match but I can't get the standard deviation in the denominator to
match. I'm not sure why.
The covariance part is correct, which is why that part matched. But the claim about the correlation matrix of $\sqrt pX+\sqrt{1-p}Y$ is incorrect, which is why the algebra doesn't work out. I will demonstrate this with simulation from the bivariate normal distribution using $p=\frac{1}{2}$. The relevant intuition in my choice of parameters is that the variances/covariance of $X$ swamp any variation in $Y$, so that the correlation structure of $\frac{1}{\sqrt{2}}X + \frac{1}{\sqrt{2}}Y$ is largely determined by $X$. Here's the R code for the set-up:
library(MASS) # has mvrnorm to simulate multivariate normals
n <- 1e7 # simulated sample size
mu <- c(0,0)
Sigma1 <- matrix(c(100,50,50,100), nrow=2)
Sigma2 <-matrix(c(1,0,0,1), nrow=2)
# covariance = 50, correlation = 0.5
x <- mvrnorm(n, mu, Sigma1)
# covariance = correlation = 0
y <- mvrnorm(n, mu, Sigma2)
# covariance = ?, correlation = ?
z <- sqrt(0.5)*x + sqrt(0.5)*y
For $X$ and $Y$ we find the covariance between the first and second components of the vector is pretty close to the specified population covariance. For $\frac{1}{\sqrt{2}}X + \frac{1}{\sqrt{2}}Y$ we find the covariance is indeed the weighted average ($p = \frac{1}{2}$) between the covariances for $X$ and $Y$.
> cov(x[,1], x[,2])
[1] 50.0234
> cov(y[,1], y[,2])
[1] -0.0004923819
> cov(z[,1], z[,2])
[1] 25.01153
For the correlation between first and second components, again the sample of simulated $X$ and $Y$ behaves in the expected manner, but $\frac{1}{\sqrt{2}}X + \frac{1}{\sqrt{2}}Y$ largely reflects the behaviour of $X$.
> cor(x[,1], x[,2])
[1] 0.5001478
> cor(y[,1], y[,2])
[1] -0.0004925274
> cor(z[,1], z[,2])
[1] 0.4951888
Let's visualise the problem by looking at the relationship between the first and second components for our samples of each of the three random vectors. To get it to plot in a reasonable time I resampled with n <- 1e5 first but this didn't substantially alter the results.
require(ggplot2)
df <- data.frame(variable=rep(c("X", "Z", "Y"), each=n),
first=c(x[,1], z[,1], y[,1]),
second=c(x[,2], z[,2], y[,2]))
g <- ggplot(df, aes(x=first, y=second, colour=variable)) +
geom_point(alpha = 0.1) + coord_fixed() + theme_bw() +
guides(colour = guide_legend(override.aes = list(alpha = 1))) +
xlab("First component") + ylab("Second component")
print(g)
g <- g + guides(colour=FALSE) + facet_grid(. ~ variable)
print(g)
We can see that the components of $Y$ are not correlated (the points from the sample of $Y$ form a roughly circular blob) while the components of $X$ are (its points form a roughly elliptical blob with its axes not aligned with horizontal and vertical). Their weighted average $Z$ resembles a scaled-down version of $X$ because the contribution from $Y$ is so small. This produces the required averaging of covariances, but leaves the correlation almost unchanged from that of $X$.
So if the original poster wants to make headway with their $\sqrt pX+\sqrt{1-p}Y$ approach (which is rather neat) then it's important to prevent such "swamping". This is feasible, since the covariance matrix does behave in the manner they desire, and a correlation matrix is just the covariance matrix for a vector whose components each have unit variance. So rather than setting $X$ as any random vector with correlation $R$, and $Y$ as any random vector independent of $X$ and with correlation $Q$, try being stricter by demanding that $Q$ and $R$ are the covariance matrices. Can you make your algebra work now? Look at this simulation with $p=\frac{1}{4}$ so that $Z$ is mostly weighted towards $Y$:
> Sigma1 <- matrix(c(1,0.4,0.4,1), nrow=2) # variance = 1 for each component
> Sigma2 <-matrix(c(1,0.8,0.8,1), nrow=2)
>
> # covariance = correlation = 0.4
> x <- mvrnorm(n, mu, Sigma1)
>
> # covariance = correlation = 0.8
> y <- mvrnorm(n, mu, Sigma2)
>
> # p = 1/4
> # covariance = ?, correlation = ?
> z <- sqrt(0.25)*x + sqrt(0.75)*y
>
> cor(x[,1], x[,2])
[1] 0.4001048
> cor(y[,1], y[,2])
[1] 0.8001744
> cor(z[,1], z[,2]) # weighted average .25*.4 + .75*.8 = 0.7
[1] 0.7000963
