The goal of this paper is to derive a formula for the finite dimensional joint characteristic function (the Fourier transform of the finite dimensional distribution) of the coupled process ${(W_{t},L_{t}^{A}):t\in \lbrack 0,\infty)}$, where $\{W_{t}:t\in \lbrack 0,\infty)}$ is a $d$-dimensional Brownian motion and $\{L_{t}^{A}:t\in \lbrack 0,\infty)}$ is the generalized $d$-dimensional L$\acute{e}$vy's stochastic area process associated to a $d\times d$ matrix $A.$ Here $A$ need not be skew-symmetric, and in our computation we allow $A$ to vary. The problem finally reduces to the solution of a recursive system of symmetric matrix Riccati equations and a system of independent first order linear matrix ODEs. As an example, the two dimensional Lévy's stochastic area process is studied in detail.