I have a functional response $Y(t)$ which I regress on a set of scalar explanatory variables $X_1$ ($=1$, i.e., $X_1$ is a constant term), $X_2$, and $X_3$. The equation I am considering is thus $$Y(t) = \alpha(t) + \beta(t)X_2 + \gamma(t)X_3 + \varepsilon(t).$$

I am using the R package fda. I know how to run the regression:

set.seed(313) # seed

t = 100 # time points
n = 50 # sample size

# generate some mock data
Y = sweep(matrix(rnorm(n * t, 0, 1), ncol = 50), 1, sin(1:100)) 
X2 = sample(c(0, 1), size = n, replace = TRUE)
X3 = rnorm(n)

# generate a basis for Y
Yfd = fda::smooth.basis(1:t, t(Y), fda::create.bspline.basis(c(1, 100), 2, 1))$fd

# run the regression
ffit = fda::fRegress(Yfd ~ X2 + X3)

Using plot(fda::predict.fRegress(ffit)) gives me plots of the fitted values. However, I want to to compute the predictions by hand. How to do that? The documentation is, unfortunately, not very helpful.

My idea was to compute $$\boldsymbol\theta(t) = (\boldsymbol X'\boldsymbol X)^{-1}\boldsymbol X'\boldsymbol y(t)$$ for each $t$, where $$\boldsymbol\theta(t) = \begin{pmatrix}\alpha(t)\\ \beta(t) \\ \gamma(t) \end{pmatrix}, \qquad\boldsymbol X = \begin{pmatrix} 1 & X_{21} & X_{31} \\ 1 & X_{22} & X_{32} \\ \vdots & \vdots & \vdots \\ 1 & X_{2n} & X_{3n} \end{pmatrix}\qquad\text{and}\qquad\boldsymbol y(t) = \begin{pmatrix} Y_1(t) \\ Y_2(t) \\ \vdots \\ Y_n(t) \end{pmatrix}.$$ Is this the correct way to do it? I don't think so, because I have never used the basis specified for the functional response...