Kohn is easily one of my favorite humans of all time, and he was a role model to whom I looked up in great admiration for most of my academic life; in fact before this site was created, I proposed that we name it after him.
However I completely disagree with the sentence that you have quoted. Keep in mind that even though the Nobel Lecture was in 1999, Kohn was born in 1923, so I was not alive for much of his life, and I don't know what possible connotations might have surrounded the word "legitimate" back in those days; but certainly the way we use the word "legitimate" nowadays, and every dictionary definition of legitimate I've seen, would indicate that he might have been speaking in hyperbole.
Let me address now, some of the specific matters in your question:
We have accurate wavefunctions for systems with far more than ${\small N=10^3}$
In this answer I recently mentioned that CCSD(T) with local orbital methods have calculated wavefunctions for systems with up to 1023 atoms; in this case it was a lipid transfer protein (PDB: 1N89) for which I'd estimate the number of electrons is about 10,000. Kohn may have written $N\simeq 10^3$ instead of $N = 10^3$, but the order of magnitude turned out to be wrong 20 years after that quote. Surely the order of magnitude will also increase again.
Accuracy
Energy differences with CCSD(T) or even LNO-CCSD(T) in the 4-zeta basis set that they used in the above example, are likely to be accurate to within 1.5 kcal/mol for a molecule like this, whereas DFT is unlikely to give you anything with an error less than 4 kcal/mol unless you use hybrids (which by definition use wavefunctions). A good energy difference doesn't necessarily mean an accurate wavefunction, but
coupled cluster wavefunctions are not bad at all (otherwise you wouldn't be able to calculate accurate properties like polarizabilities).
Storage
The wavefunction in the above case is stored via cluster amplitudes, for which we can store billions (in fact trillions in the largest cases)
of them. Since cluster amplitudes appear in the argument of an exponential, we actually get CI coefficients for 100% of the non-zero determinants. So storage of a big wavefunction is no problem when you use a compact representation. Sure there will be a point at which classical computers can no longer accurately store quantum wavefunctions, but there will also be a point at which electron densities can't be stored either, so in that sense why not call the density "illegitimate" too? Furthermore, not being able to "store" the wavefunction is only a problem if using a classical storage device, whereas if you use qubits instead of bits the statement no longer has legs.
Quantum Monte Carlo (QMC)
FCIQMC, VMC, DMC, AFQMC, and similar methods, are all wavefunction methods. They are wavefunctions methods, whereas you wrote that they "could be called wave function methods, in that...". They do in fact manage to represent wavefunctions in huge Hilbert spaces, often by taking advantage of sparsity of the wavefunction, but using a compact representation as in the case of coupled cluster, means that you don't even need a "sparse" representation or a "stochastic sampling", you can represent the whole wavefunction by storing only the argument of the exponential. I am not trying to take shine away from QMC; I say all this as a contributor to a major FCIQMC code, and having used FCIQMC in papers with the inventor of FCIQMC and also separately on my own and even on this paper I put on arXiv only two days ago. FCIQMC has its place as one of the best methods for ultra-high accuracy in large multi-reference systems, but it is not needed in the proof that Kohn's statement is wrong: coupled cluster can be very accurate on even bigger systems if they don't have too much multi-reference character.
Number of electrons is actually a red herring
The problem with wavefunctions when using classical computers has less to do with the number of electrons and more to do with its structure:
- Bosonic wavefunctions don't suffer from the fermionic sign-problem, so you can represent them with a Hartree product, and might be able to store a larger wavefunction than you can store a density! Kohn's statement is about electronic wavefunctions, but this bullet serves as a reminder that Kohn's argument is not so much about fundamental physics and ontology, than it is about "computability" (you probably already agreed with this, so this bullet is more for other people).
- Fermionic wavefunctions involving only one determinant (which can still be accurate for a very single-reference system) are very simple: instead of the $\binom{M}{N}$ type scaling for the number of determinants in a full CI expansion, you only have one term.
- Fermionic wavefunctions involving many determinants but only static correlation, can be represented by matrix product states and calculated using the polynomially scaling DMRG.
- Fermionic wavefunctions involving many determinants but only dynamic correlation, can be represented by coupled cluster ansatze which are also polynomially scaling.
- Fermionic wavefunctions involving infinitely many determinants can also be represented by a compact representation.
The problem is more about how many digits you need for each of the CI coefficients. Then you quickly see that the computational complexity of $2^N$ vs $N^3$ is irrelevant and what really matters is something more subtle, which is: how complicated is the wavefunction, not how many electrons are there.
50 electrons in a CAS(50,50) is currently an absolutely brutal calculation but 10,000 electrons in a CCSD(10000,44000) was done easily in the paper listed above. So there are cases where $N=50$ electrons is harder than $N=10^4$, and in those highly multi-reference cases, good luck getting an accurate energy with a single-reference method like DFT!
