We present a methodology for estimating causal functional linear models using orthonormal tensor product expansions. More precisely, we estimate the functional parameters α and β that appear in the causal functional linear regression model:where supp β ⊂ T, and T is the closed triangular region whose vertexes are (a, a), (b, a) and (b, b). We assume we have an independent sample {(Y k , X k ) : 1 ≤ k ≤ N, k ∈ N} of observations where the X k 's are functional covariates, the Y k 's are time order preserving functional responses and E k , 1 ≤ k ≤ N, is i.i.d. zero mean functional noise.