Objective
Let's say that \$I_{_\text{C}}\$ is the DC operating point value you want to evaluate for its stability against certain changes. And let's say that temperature is one of these and that device variations are another of these.
Detailing the objective
The first thing you need to do is narrow down what temperature affects and, if possible, how it affects those things. Temperature will impact both \$\beta\$ and also \$V_{_\text{BE}}\$.
There are complex equations, depending upon something called the Boltzmann factor, that impact the saturation current model parameter, which itself impacts \$V_{_\text{BE}}\$. The range seems to be about from \$-1.8\:\text{m}\:\frac{\text{V}}{^\circ\text{C}}\$ to about \$-2.4\:\text{m}\:\frac{\text{V}}{^\circ\text{C}}\$ But in general, we can say that \$V_{_\text{BE}}\$ varies by about \$-2\:\text{m}\:\frac{\text{V}}{^\circ \text{C}}\$ near room temps. This works out to about a \$- 1\% \$ change for each \$+ 1\% \$ change in temperature.
\$\beta\$ will be a different thing. It increases with increases in temperature. But let's just assume for now that \$\beta\$ changes by \$+ 1\% \$ for each \$+ 1\% \$ change in temperature.
The above doesn't account for changes in the Early Effect and a variety of other impacts. But we can assume those away, for now.
Capturing the boundaries
So, for a \$\pm 30\:^\circ\text{C}\$ temperature range we can expect to see \$V_{_\text{BE}}\$ varies about \$\pm 10\% \$, but inversely so, and \$\beta\$ varies also about \$\pm 10\% \$, but in the same direction.
As regards device variations, let's say that \$V_{_\text{BE}}\$ varies about \$\pm 10\% \$ and that \$\beta\$ varies about \$\pm 50\% \$.
Let's assume all these elements are orthogonal to each other and apply RMS to them to get a \$V_{_\text{BE}}\$ variation of about \$\pm 15\% \$ and a \$\beta\$ variation of about \$\pm 50\% \$.
This gives us some magnitudes to work with for the circuit.
Quantitative circuit analysis: \$I_{_\text{C}}\$
Oh. Almost forgot. What about the circuit? I guess we need some formulas!
If you analyze your circuit using KCL and KVL you will find that:
$$I_{_\text{C}}=\beta\cdot\frac{V_{_\text{CC}}-V_{_\text{BE}}}{R_{_\text{B}}+\left(\beta+1\right)R_{_\text{C}}}$$
To see how, read the following SymPy/Python code:
var( 'rb rc vb vc ib ic beta ie vbe vcc' )
eq1 = Eq( vc/rc + vc/rb + ic, vb/rb + vcc/rc ) # KCL collector node
eq2 = Eq( vb/rb + ib, vc/rb ) # KCL base node
eq3 = Eq( vcc - rc*(ic+ib) - rb*ib - vbe, 0 ) # KVL through resistors to ground
eq4 = Eq( ic, beta*ib ) # BJT
eq5 = Eq( ie, ic + ib ) # BJT
ans = solve( [ eq1, eq2, eq3, eq4, eq5 ], [ vc, vb, ic, ib, ie ] )
ans[ic]
-beta*(vbe - vcc)/(beta*rc + rb + rc)
pprint( ans[ic] )
-β⋅(vbe - vcc)
───────────────
β⋅rc + rb + rc
Sensitivity
A sensitivity equation is used to relate the %-change of one thing to the %-change of another thing. This looks like: \$\frac{\frac{\text{d}\, y}{y}}{\frac{\text{d}\, x}{x}}=\frac{\text{d}\, y}{\text{d}\, x}\cdot\frac{x}{y}\$, which gives you the %-change of y as compared to the %-change of x.
In this case, I find:
$$\begin{align*}
\left|\%\,I_{_\text{C}}\right|&=\%\,V_{_\text{BE}}\cdot\frac{V_{_\text{BE}}}{V_{_\text{CC}}-V_{_\text{BE}}}
\\\\
\left|\%\,I_{_\text{C}}\right|&=\%\,\beta\cdot\frac{R_{_\text{B}}+R_{_\text{C}}}{R_{_\text{B}}+\left(\beta+1\right)R_{_\text{C}}}
\end{align*}$$
The SymPy code looks like this:
simplify( derivative( ans[ic], vbe ) * vbe / ans[ic] )
vbe/(vbe - vcc)
simplify( derivative( ans[ic], beta ) * beta / ans[ic] )
(rb + rc)/(beta*rc + rb + rc)
Getting more specific
So, let's assume \$I_{_\text{C}}=2\:\text{mA}\$ and \$V_{_\text{CC}}=10\:\text{V}\$. Assuming also a design value of \$\beta=200\$ and \$V_{_\text{BE}}=700\:\text{mV}\$ and a quiescent collector operating point near \$3.4\:\text{V}\$, we might then specify \$R_{_\text{C}}=3.3\:\text{k}\Omega\$ and \$R_{_\text{B}}=270\:\text{k}\Omega\$. (Recomputing, we'd then find \$I_{_\text{C}}=1.99\:\text{mA}\$. Close enough.)
From these values and the above %-change estimates, we'd say that the collector current might change by about \$\pm\,14.6\%\$ over part and temperature variations of \$\beta\$ and by about \$\pm\,1.1\%\$ over part and temperature variations of \$V_{_\text{BE}}\$.
The SymPy code looks like:
abs( simplify( derivative( ans[ic], beta ) * beta / ans[ic] ).subs(
{ beta:200, rb:270e3, rc:3.3e3 }
) * 50 )
14.6415943426551
abs( simplify( derivative( ans[ic], vbe ) * vbe / ans[ic] ).subs(
{ beta:200, rb:270e3, rc:3.3e3, vcc:10, vbe:.7 }
) * 15 )
1.12903225806452
If the argument so far has been soundly reasoned, then it's now obvious that \$\beta\$ variation swamps out \$V_{_\text{BE}}\$ variation and you should figure that there will be as much as \$\pm\,15\:\%\$ variation in the operating point over temperature and device variations.
(You could remove the orthogonality assumption I made and argue for still worse. But this gives an idea about an approach to quantify variations you care about.)
This particular circuit topology
This circuit topology is never used in isolation where a specific voltage gain is required, because the voltage gain is highly dependent upon the collector current and because it is also highly dependent upon the signal, itself. So this is used to get the highest gain possible from a BJT, with the expectation of very high distortion. It is therefore always wrapped with negative feedback by placing it within a larger circuit context to remove distortion and set the desired voltage gain.
The analysis so far tells us something useful about this topology.
To get the highest possible voltage gain from this circuit, \$R_{_\text{C}}\$ should be as large as possible without forcing the BJT into deep saturation. But we must allow for collector current variation of about \$15\,\%\$. We might have been tempted to set \$R_{_\text{C}}=3.9\:\text{k}\Omega\$ to get better voltage gain. But then we'd find a quiescent collector voltage that could be under one volt. And that's without any signal applied! So, with signal, we'd almost certainly be driving the BJT well into deep saturation. That's why this kind of analysis can be vital in deciding how much to compromise on the voltage gain in order to deal with circuit variations, even with NFB applied.
Engineering logbooks
As an aside, what you read above is something that should be placed into the engineering logbook you keep as a circuit designer. (A must-do, if you are a professional, and should-do, otherwise.) These kinds of details are important to record and to be able to refer back to, should circumstances require it, in order to show why you made the choices you made at the time. For professionals, it can save your life if things ever go to court. For professionals and hobbyists, it can save a lot of your time later on when you need to revisit a design.
In short, keep logbooks for the rest of your life!