Expected number of cases, standardizing with age & sex of the whole population Starting with a toy example data (http://www.stata-press.com/data/r14/hbp.dta)
use http://www.stata-press.com/data/r14/hbp, clear

drop if mi(sex)
drop if mi(age_group)
sort city sex age_group year, stable
by city sex age_group: egen obs = total(hbp)
by city sex age_group:  gen pop = _N
by city sex age_group:  keep if _n == 1
drop year race hbp

replace city = 4 if city == 5 // used by the the loop later, so sequential better

bysort city (sex age): egen obs_city = total(obs)
bysort city (sex age): egen pop_city = total(pop)
gen crude = obs_city / pop_city

We have age & sex specific counts of events obs and population denominator pop in 5 different cities. Crude rate in a city is given by crude variable simply dividing observations by population within city and ignoring age structure.  Now, I'd like to use information from all cities to create age & sex specific rates and using that - calculated expected number of events given the study population. 
I'm trying to achieve that by using city specific standardized rates:
dstdize obs pop age sex, by(city) print saving( "temp", replace )

mat b = r(adj)
matrix list b

gen adjusted = .
gen exp_city_ds = .

forv city = 1/4 {
    local bw = b[1, `city']
    display `bw'
    replace adjusted = `bw' if city == `city'
    replace exp_city_ds = pop_city * `bw' if city == `city'
}

So far so good, but when I look at results I'd obtain from simple poisson model - the numbers do not add up:
poisson obs i.sex i.age, exp(pop) irr 
predict exp_p, n

bysort city (sex age): egen exp_city_poi = total(exp_p)
drop exp_p

gen dif0 = obs_city - exp_city_ds   
gen dif1 = obs_city - exp_city_poi  
gen dif2 = exp_city_ds - exp_city_poi
su dif? if sex == 1 & age == 2 // first obs by city

saveold "temp_old.dta", v(12) replace // data for R

I also tried to calculate from the same dummy dataset observed number of case using expected function of SpatialEpi package:
library(foreign)
library(SpatialEpi)

data <- read.dta("temp_old.dta") #, convert.factors = FALSE)

expected(data$pop, data$obs, 8)
seq(1, 4, 1)

city <- cbind(seq(1, 4, 1), expected(data$pop, data$obs, 8))
colnames(city) <- c("city", "E")

data <- merge(data, city, by="city")

data$poi_fit <- fitted(glm(obs ~ 1 + offset(log(pop)), data = data, family = "poisson"))

Again - I get different numbers, very close to Poisson fit (both Stata & R) but not the same.
Now, my question: is the reason of such difference lying in my calculations/methods? Or is it expected?
 A: You are applying Poisson model to the dataset which does not behave this way -- the assumption of independence in time does not hold.
For the dataset in consideration. We look at the prevalence of hypertension for different cities and age/sex groups.
Try the following R code:
# As example read high blood pressure data
library(foreign)
d <- read.dta("hbp.dta")
# Reformat data
d$city <- as.factor(d$city)
d$age_group <- factor(d$age_group) # Remove misleading unused age groups
d$hbp <- as.numeric(d$hbp) - 1
# Fit logistic regression
fm <- glm(hbp ~ city + age_group*sex, data=d)
# Basic model selection
bm <- step(fm)
summary(bm)

You will find that the risk age group turns out to be (30-34) which are, indeed, the oldest subjects in the data. Being a woman is slight protective factor and the city3 has the highest prevalence. Therefore to test the Poisson behavior let us subset to, say, city1, men.
ind1 <- which(d$city == "1")
ind2 <- which(d$sex == "Male")
d2 <- d[c(ind1,ind2),]
ag <- aggregate(hbp ~ age_group, d2, sum)

The result is:
age_group hbp
15 - 19   5
20 - 24  12
25 - 29  19
30 - 34  24

Now fit Poisson:
x <- (1:4)
m1 <- glm(ag$hbp ~ x, family="poisson")

The model has $R^2 = 0.96$ so it is good.
