You must use the standard deviation of the prediction errors, rather than the standard deviation of the residuals. The former varies for each forecast period.
For example, in an AR(1) model
$$y_t=\phi y_{t-1} + \varepsilon_t \,, \quad \varepsilon_t\sim NID(0, \sigma^2) \,, \quad t=1,\dots,T\,,$$
the variance of the prediction errors can be obtained as follows.
The h-steps ahead prediction errors are defined as $e(h) = y_{T+h} - E(y_{T+h})$ (the value determined by the model minus its expectation):
$$
\begin{align}
e(1) &= y_{T+1} - E(y_{T+1}) = (\phi y_T+\varepsilon_{T+1}) - \phi y_T = \varepsilon_{T+1} \\
e(2) &= y_{T+2} - E(y_{T+2}) = (\phi \underbrace{y_{T+1}}_{\phi y_T+\varepsilon_{T+1}}+\varepsilon_{T+2}) - E(y_{T+2}) \\
&= \underbrace{(\phi^2y_T+\phi \varepsilon_{T+1}+\varepsilon_{T+2})}_{y_{T+2}} - \phi^2y_T) \\
&= \phi \varepsilon_{T+1} + \varepsilon_{T+2} \\
e(3) &= \underbrace{\phi^3y_T+\phi^2\varepsilon_{T+1}+\phi \varepsilon_{T+1} + \varepsilon_{T+3}}_{y_{T+3}} - E(y_{T+3}) \\
&= \phi^3y_T+\phi^2\varepsilon_{T+1}+\phi \varepsilon_{T+1} + \varepsilon_{T+3} - \phi^3y_T = \phi^2\varepsilon_{T+1}+\phi \varepsilon_{T+2} + \varepsilon_{T+3} \\
\cdots & = \cdots \\
e(h) &= \phi^{h-1}\varepsilon_{T+1}+\dots+\phi \varepsilon_{T+h-1} + \varepsilon_{T+h}
\end{align}
$$
The variances are straightforward to obtain (just note that $Var(\phi^k \varepsilon_t)=\phi^{2k}\sigma^2$):
$$
\begin{align}
Var(e(1)) &= \sigma^2 \\
Var(e(2)) &= \sigma^2(1+\phi^2) \\
Var(e(3)) &= \sigma^2(1+\phi^2+\phi^4) \\
\dots &= \dots
\end{align}
$$
The expressions of the variances vary with the order of the model.
A common way to compute the variances of the prediction errors in a general ARMA(p,q) is by means of the Kalman filter. Here is some code for illustration.
These are your data:
x <- ts(c(319.32, 320.36, 320.82, 322.06, 322.17,321.95,321.20,318.81,317.82, 317.37, 318.93, 319.09, 319.94, 320.98, 321.81, 323.03,323.36, 323.11, 321.65, 319.64, 317.86, 317.25, 319.06, 320.26, 321.65, 321.81,322.36, 323.67, 324.17, 323.39, 321.93, 320.29, 318.58, 318.60,319.98, 321.25, 321.88, 322.47, 323.17, 324.23, 324.88, 324.75, 323.47, 321.34,319.56, 319.45, 320.45, 321.92, 323.40, 324.21, 325.33, 326.31,327.01, 326.24, 325.37, 323.12, 321.85, 321.31, 322.31, 323.72, 324.60, 325.57,326.55, 327.80, 327.80, 327.54, 326.28, 324.63, 323.12, 323.11, 323.99, 325.09, 326.12, 326.61, 327.16, 327.92, 329.14, 328.80, 327.52,325.62, 323.61, 323.80, 325.10, 326.25, 326.93, 327.83, 327.95,329.91, 330.22, 329.25, 328.11, 326.39, 324.97, 325.32, 326.54, 327.71, 328.73, 329.69, 330.47, 331.69, 332.65, 332.24, 331.03, 329.36, 327.60, 327.29, 328.28, 328.79, 329.45, 330.89, 331.63, 332.85, 333.28, 332.47, 331.34, 329.53, 327.57, 327.57, 328.53, 329.69, 330.45, 330.97, 331.64, 332.87, 333.61, 333.55, 331.90, 330.05, 328.58, 328.31, 329.41, 330.63, 331.63, 332.46, 333.36, 334.45, 334.82, 334.32, 333.05, 330.87, 329.24, 328.87, 330.18, 331.50, 332.81, 333.23, 334.55, 335.82, 336.44, 335.99, 334.65, 332.41, 331.32, 330.73, 332.05, 333.53, 334.66, 335.07, 336.33, 337.39, 337.65, 337.57, 336.25, 334.39, 332.44, 332.25, 333.59, 334.76, 335.89, 336.44, 337.63, 338.54, 339.06, 338.95, 337.41, 335.71, 333.68, 333.69, 335.05, 336.53, 337.81, 338.16, 339.88, 340.57, 341.19, 340.87, 339.25, 337.19, 335.49, 336.63, 337.74, 338.36), frequency=12, start=c(1965,1))
These are the results from forecast that you want to reproduce:
require("forecast")
fit <- Arima(x, order = c(1,1,0), seasonal = list(order =c(0,1,1), period =12))
forecast(fit, h=4)
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
# Jan 1981 339.5735 339.1483 339.9987 338.9232 340.2238
# Feb 1981 340.1759 339.6493 340.7025 339.3706 340.9813
# Mar 1981 341.2141 340.5865 341.8417 340.2543 342.1739
# Apr 1981 342.3020 341.5914 343.0125 341.2153 343.3886
The state space representation of the fitted model is available in fit$model, which is the input required to run the Kalman filter.
kf <- KalmanForecast(n.ahead = 4, fit$model)
pred.vars <- kf$var * fit$sigma2
kf$pred - 1.96 * sqrt(pred.vars)
# [1] 338.9232 339.3706 340.2543 341.2153
kf$pred + 1.96 * sqrt(pred.vars)
# [1] 340.2238 340.9813 342.1740 343.3886
The variances of prediction errors returned by the Kalman filter are scaled by the estimated variance of the residuals, $\sigma^2$. Then, the confidence intervals are obtained as you did (you can see that these values match the columns Lo 95 and Hi 95 obtained above).
Edit
This may depart from your question, but as mentioned by @IrishStat, if the assumption of Gaussian innovations is not plausible, it is worth considering some alternative methods. For example, the option bootstrap=TRUE computes bootstrapped confidence intervals. After a quick view to the documentation, I think this option resamples the residuals and generates an ensemble of replicates based on the parameters of the fitted model. Although the distribution of the residuals of your fitted model depart from the Gaussian distribution (due to excess of kurtosis as suggested for example by the histogram), graphically there are no major differences in the confidence bands. Here are the numbers:
set.seed(1)
forecast(fit, h=4, bootstrap=TRUE)
# Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
# Jan 1981 339.5735 339.2259 339.9744 338.9472 340.2538
# Feb 1981 340.1759 339.7288 340.6966 339.4446 341.0239
# Mar 1981 341.2141 340.6837 341.8561 340.3736 342.2230
# Apr 1981 342.3020 341.6879 343.0394 341.3340 343.4285
Aside note: I think your question is more related to theoretical issues
rather than on practical issues and your are using this model and data
mainly for illustration. However, as preparing this answer I had a closer look to your data, I would suggest an ARIMA(0,1,1)(0,1,1) model with a temporary change at observation 190 (october 1980). I don't claim it is the best model or that it will perform better for forecasting, but according to the diagnostic of the residuals, it is better than the ARIMA(1,1,0)(0,1,1).
