There was a slight error in part (ii) of my question. I have made the change: $M(\beta)$ should be been $M^\prime(\beta)$, it's derivative.
Here is my own approach: please provide your comments.
For all $(y,x)\in\{0,1\}\times \mathcal{X}$, it is easy to check that the map, $\beta\mapsto y\log{\psi(x\beta)}+(1-y)\log{\{1-\psi(x\beta)\}}$, is strictly concave, where $\mathcal{X}$ is the sample space for $X$. I am assuming $\beta\in\mathbb{R}$. Let, $P_0$ denote the joint d.f. of $(Y,X)$ under $\beta_0$. As $Q(x)$ is non-degenerate, $P_0$ will be non-degenerate and hence,
$$
\beta\mapsto M(\beta) = \int \left[y\log{\psi(x\beta)}+(1-y)\log{\{1-\psi(x\beta)\}}\right]~dP_0(y,x),
$$
will also be strictly concave(?). Now, $M(\beta)\leq 0$, for all $\beta$, since $\psi(\cdot)\in (0,1)$.
My doubt: Is it possible that a strictly concave map defined on the real line can take the value $-\infty$, at some point $u\in \mathbb{R}$? I need some help on this.
Now I assume: $M(\beta)>{}-\infty$, for all $\beta$. So, at any $\beta\in\mathbb{R}$, and any $a_n\rightarrow 0$, we have,
\begin{align*}
&\lim_{a_n\rightarrow 0} \frac{M(\beta + a_n) - M(\beta)}{a_n} \\
&=\lim_{a_n\rightarrow 0} \int y\frac{[\log{\psi(x\beta + x a_n)}-\log{\psi(x\beta)}]}{a_n}~dP_0\\
&\ {}+\lim_{a_n\rightarrow 0} \int (1-y)\frac{\big[\log{\{1-\psi(x\beta + x a_n)\}}-\log{\{1-\psi(x\beta)\}}\big]}{a_n}~dP_0
\end{align*}
Consider the first term and write, $g_1(u) = \log{\psi(u)}$. Using Mean value theorem, at each fixed $(y,x)$, we have
$$
g_1(x\beta + xa_n) - g_1(x\beta) = (xa_n)\cdot g^\prime_1(x\beta + \eta_{x,n} (xa_n)),\quad\text{for some $\eta_{x,n}\in (0,1)$.}
$$
Here, $g^\prime_1(u) = 1- \psi(u)$. Hence, the first term on the rhs becomes,
\begin{align*}
\lim_{a_n\rightarrow 0} \int yx\big\{1 - \psi(x\beta + \eta_{x,n}(xa_n))\big\}~dP_0 = \lim_{n\rightarrow\infty}\int g_n(y,x)~dP_0,
\end{align*}
where, $g_n(y,x)$ is the function inside the integral. Note, $|g_n(y,x)|\leq 2|x|$, which is $Q(\cdot)$-integrable, by assumption. Also, for each $(y,x)$, $g_n(y,x)\rightarrow yx\{1-\psi(x\beta)\}$. Now, using Dominated convergence theorem, the first term on the rhs converges to $\int x y\{1-\psi(x\beta)\}~dP_0$. Similarly, the second term converges to $\int x (1-y)\{-\psi(x\beta)\}~dP_0$. Eventually, we obtain $M^\prime(\beta) = \mathbf{E}\big(Z(\beta:Y,X)\big)$.
My doubt: The way I have used Mean value theorem inside integral, at each $(y,x)$. Is it rigorous, or have I missed anything.
I will highly appreaciate your comments.