Until the recent advent of $Gaia$ Data Release 2 (DR2) and deep multi-object spectroscopy, it has been difficult to obtain 6-D phase space information for large numbers of stars beyond 4 kpc, in particular towards the Galactic centre, where dust and crowding effects are significant. In this study we combine line-of-sight velocities from the Abundances and Radial velocity Galactic Origins Survey (ARGOS) spectroscopic survey with proper motions from $Gaia$ DR2, to obtain a sample of $sim$ 7,000 red clump stars with 3-D velocities. We perform a large scale stellar kinematics study of the Milky Way (MW) bulge to characterize the bulge velocity ellipsoids. We measure the tilt $l_{v}$ of the major-axis of the velocity ellipsoid in the radial-longitudinal velocity plane in 20 fields across the bulge. The tilt or vertex deviation, is characteristic of non-axisymmetric systems and a significant tilt is a robust indicator of non-axisymmetry or bar presence. We compare the observations to the predicted kinematics of an N-body boxy-bulge model formed from dynamical instabilities. In the model, the $l_{v}$ values are strongly correlated with the angle ($alpha$) between the bulge major-axis and the Sun-Galactic centre line-of-sight. We use a maximum likelihood method to obtain an independent measurement of $alpha$, from bulge stellar kinematics alone. The most likely value of $alpha$ given our model is $alpha = (29 pm 3)^{circ}$. In the Baades window, the metal-rich stars display a larger vertex deviation ($l_{v} = -40^{circ}$) than the metal-poor stars ($l_{v} = 10^{circ}$) but we do not detect significant $l_{v}-$metallicity trends in the other fields.